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severely  Indented  nosetip  shapes  under  investigation,  routine  calculation  fails 
to  give  a  solution  and  a  special  calculation  procedure  has  been  developed  in 
order  to  obtain  reasonable  solutions.  For  viscous  flows,  extensive  calculations 
have  been  performed  for  both  smooth  and  indented  nosetip  shapes,  however, 
comparisons  of  calculated  results  and  measured  data  are  not  satisfactory.  For 
smooth  nosetip  shapes,  although  the  surface  pressure  and  shock  location  agree 
well  but  the  temperature  field  is  poor.  For  indented  nosetip  shapes, 
difficulties  to  simulate  flows  with  large  separation  bubble  and  sharp  corner 
as  well  as  turbulent  flows  are  described. 

Because  of  the  poor  results  obtained  for  viscous  flow  calculations,  a 
reanalysis  of  the  complete  calculation  procedure  for  solving  the  full  Navier- 
Stokes  equations  has  been  carried  out.  The  original  code  was  then  modified  by 
rewritting  all  the  viscous  subroutines  according  to  the  new  analysis  for  thin- 
layer  theory  without  turbulence  model.  The  modified  nosetip  code  has  been 
applied  to  compute  laminar  flow  over  a  hemisphere-cylinder  and  a  sphere-cone. 

The  new  results  do  not  show  that  temperature  is  a  mesh  dependent  variable  as 
reported  by  Kutler  et  al  and  the  temperature  field  compare  well  with  the 
available  solution  for  full  Navier-Stokes  equations  and  boundary  layer 
calculation. 

Finally,  a  coupling  of  the  nosetip  code  with  an  existing  NSWC  supersonic 
marching  code  for  inviscld  afterbody  solution  has  been  carried  out.  Results 
for  surface  pressure  are  presented  for  a  very  blunt  nose  cone  and  an  indented 
nosetip  cone  and  the  agreement  is  good. 
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FOREWORD 


A  finite  difference  computer  program  has  been  developed  at  NSWC  to  calculate 
the  inviscid  and  viscous  flow  fields  for  blunt  or  indented  nosetips  for  reentry 
vehicles  at  hypersonic  flight.  This  report  describes  the  development  of  the 
computer  program,  illustration  of  numerical  examples  and  the  operating  instruc¬ 
tions. 

This  work  was  supported  by  the  Reentry  Aerodynamic  Program  of  NSWC  at  White 
Oak  and  monitored  by  Drs.  A.  M.  Morrison  and  W.  C.  Lyons. 

The  author  would  like  to  thank  Dr.  Paul  Kutler  of  NASA  Ames  Research  Center 
for  providing  a  copy  of  the  -computer  code  and  the  instructions  for  its  operation. 
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INTRODUCTION 

Because  of  ablation,  the  nosetip  of  a  spherical  body  undergoes  continuous 
change  during  re-entry.  The  shape  of  the  nosetip  has  a  great  influence  on  the 
flowfield  over  the  entire  body,  i.e.,  the  nose  region  and  thus  the  afterbody.  In 
order  to  predict  the  flowfield  about  Indented  nosetip  shapes  that  are  likely  to 
occur  during  the  reentry,  considerable  effort  has  been  expended  in  the  numerical 
simulation  of  hypersonic  flow  over  indented  nosetips.l“5 


Among  the  many  numerical  schemes  intended  for  indented  nosetip  calculation,  an 
attractive  one  seems  to  be  due  to  Kulter  et  al , *  who  solve  the  unsteady  Navier- 
Stokes  equations  with  the  thin-layer  approximation  for  nosetip  of  arbitrary  shapes 
at  zero  incidence  using  the  implicit  factored  numerical  algorithm  of  Uarming  and 
Beam.  The  steady  solution  is  obtained  asymptotically  in  time  and  both  viscous 
and  inviscld  flowfields  can  be  computed  using  the  same  computer  program. 


The  obvious  reason  for  choosing  implicit  scheme  is  that  it  is  more  efficient 
for  viscous  flow  calculation  and  one  expects  flow  separation  to  play  an  important 
role  in  nosetip  flowfield  simulation.  A  research  code  was  then  obtained  from 
Dr.  Kutler  of  NASA  Ames  Research  Center.  First,  the  code  was  applied  to  compute 
inviscld  flow  over  sphere  (or  sphere-cone),  results  compared  well  with  experiments. 
When  the  similar  calculations  were  performed  for  a  series  of  four  indented 
nosetip  shapes  reported  in  Ref.  7  and  8,  surmountable  difficulties  were  encountered 
during  the  course  of  computation  because  of  the  presence  of  small  radius  expandion 
corner  and  concave  compression  turn  in  these  nosetips.  It  was  later  found  that  a 
special  calculation  procedure  is  required  in  order  to  obtain  reasonable  solutions. 
Example  to  demonstrate  this  special  calculation  procedure  will  be  given  and 
the  code  was  modified  accordingly  to  do  the  job.  All  inviscld  results  for  indented 
nosetips  are  described  in  Section  5. 
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Later,  viscous  laminar  flow  over  a  sphere-cone  with  isothermal  wall  was 
calculated  and  results  for  heat  transfer  rate  compare  reasonably  well  with 
measured  data  (within  experimental  error  band)  but  are  higher  than  that  predicted 
by  boundary  layer  calculation.  Then,  calculations  were  performed  for  laminar  flow 
over  two  of  the  indented  nosetip  shapes  under  investigation  using  the  same  special 
calculation  procedures  as  used  for  inviscid  calculations.  Solutions  obtained  were 
not  satisfactory  as  compared  to  the  measured  data  for  surface  pressure  and 
separation  region  (also  negative  heat  transfer  rate  were  found  in  part  of  the  body 
surface).  For  turbulent  flow  calculation  with  the  algebraic  turbulence  model  of 
Baldwin  and  Lomax,  it  was  found  not  possible  to  obtain  a  converged  solution  for 
these  indented  nosetip  shapes.  Finally,  attempt  was  made  to  repeat  the  hemisphere- 
cylinder  results  given  in  Figs.  5  to  8  of  Ref.  1  and  was  not  able  to  duplicate  the 
temperature  field  as  shown  in  Fig.  7  of  Ref.  1.  With  all  the  difficulties  in  the 
viscous  flow  calculations,  which  are  described  in  detail  in  Section  6,  it  was 
decided  that  a  reanalysis  of  the  complete  calculation  procedure  must  be  performed 
and  all  the  viscous  subroutines  must  be  rewritten.  The  analysis  is  given  in 
Sections  1  and  2.  The  boundary  conditions  given  in  Section  3  are  the  same  as 
described  in  Ref.  1,  but  the  initial  condition  has  been  modified. 


The  modified  computer  code  has  since  been  applied  to  compute  laminar  flows  over 
a  hemisphere-cylinder  with  adiabatic  wall  and  a  sphere-cone  with  isothermal  wall. 

New  results  for  the  temperature  field  of  the  hemisphere-cylinder  case  are  compared 
to  the  available  solution  given  by  Viviand  and  Ghazzi,  ^  who  solved  the  full 
Navier-Stokes  equations  and  that  given  by  Kutler  et  al.  The  dubious  statement 
that  temperature  is  a  mesh  dependent  variable  as  given  in  Ref.  1  does  not  appear 
in  the  present  results.  New  results  for  the  heat  transfer  rate  in  terms  of 
Stanton  number  for  the  sphere-cone  case  are  compared  to  the  result  obtained  with 
original  code,  the  experimental  data  and  a  boundary  layer  calculation.  Surprisingly 
good  agreement  is  shown  for  the  heat  transfer  rate  between  the  thin-layer  and  the 
boundary  layer  calculations  as  expected.  All  the  new  results  are  described  in 
Section  7.  The  entire  code  has  since  been  rewritten  for  the  convenience  of 
applying  to  indented  nosetip  shapes  described  in  this  report.  Simple  operational 
manual  and  the  computer  program  are  given  in  Appendix  D.  Due  to  the  termination 
of  support,  no  further  calculation  has  been  made  for  the  indented  nosetip  shapes 
with  the  new  code.  Also  no  investigation  of  the  behavior  of  turbulence  models 
with  the  new  code  has  been  performed. 


Since  the  nosetip  flowfield  calculation  is  to  provide  a  starting  solution  for 
the  afterbody  calculation.  Thus,  a  coupling  of  the  nosetip  code  with  an  existing 
NSWC  supersonic  marching  code  for  inviscid  calculation  has  been  carried  out. 
Examples  are  given  for  a  very  blunt  cone  at  =  9.8  and  an  indented  nosetip  cone 
at  Hd  —  5 . 
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CHAPTER  1 

GOVERNING  EQUATIONS 


The  time-dependent  compressible  Navier-Stokes  equations  in  the  cylindrical 
coordinates  (x,y,$)»  Fig.  1,  for  axisymmetric  flow  can  be  written  in  dimensionless 
conservation- law  form  for  a  perfect  gas  without  external  force  as  follows^: 


where 


and 
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Peyret,  R. ,  and  Viviand,  H.,  "Computations  of  Viscous  Compressible  Flows  Based 
on  the  Navier-Stokes  Equations,"  AGARD-AG-212,  Sep  1975. 
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A  BOW  SHOCK 
B  INTITIAL  BODY 
B'  FINAL  BODY 
C  AXIS  OF  SYMMETRY 
D  OUTFLOW  BOUNDARY 


b.  COMPUTATIONAL  PUNE 


FIGURE  1.  COORDINATE  SYSTEM 


In  Eq.  (1.1),  the  reference  quantities  used  to  nondimensionalize  the  flow 
variables  are:  the  length  L,  the  velocity  a<»//f,  the  density  Poo»  the  viscosity  y<» 
and  the  thermal  conductivity  The  time  is  nondinemsionalized  by  /7  L/aoo,  the 

total  energy  per  unit  volume  e  andpressure  p  by  poo  (=P<»ao£/Y)  and  the  viscous 
stress  term  Txx  etc.  by  a,*,  /(l/7).  For  perfect  gas,  the  equation  of  state  gives 

p  *  (Y  -  1)  [e  -  *sp(u2  +  v2)]  (1.2) 

and 

e  =  p[e  +  %  (u^  +  v^)]  (1.3a) 

e  =  Cv  T  (1.3b) 

Equations  (1.1)  and  (1.2)  provide  5  equations  for  5  unknowns,  i.e.  p,  u,  v,  e,  p 
(or  T).  It  should  be  pointed  out  that  the  momentum  equations  and  the  energy 
are  of  parabolic  type  and  the  continuity  equation  carrys  the  hyperbolic  character. 
Therefore,  the  system  of  the  N.  S.  equations  is  of  hybrid  parabolic  and  hyperbolic 
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type  [12]. 

A  mapping  between  the  physical  plane  and  the  computational  plane  (Fig.  1) 
is  accomplished  by  the  following  independent  variable  transformation: 

T  =  t,  5  *  £  (t,x,y),  T]  =  n  (t.x.y)  (1.4) 

Applying  Eq.  (1.4)  to  Eq.  (1.1),  one  obtains 

Ur  +  E5  +  Fn  +  H  -  <R5  +  S„  +  T)  d-5) 

where 


U  -  U/J 

E  =  <5t  0  +  5X  E  +  Cy  F>/J 

F  =  (nt  u  +  nx  E  +  \  F)/J 

h  -  (f  +  H)/(yJ) 

R  -  (Sx  R  +  Cy  S)/J 
s  =  (nx  R  +  Hy  S)/J 
T  *  (S  +  T)/(yJ) 

with  metrics  of  transformation  given  by 

h  -  <xn  yr  -  y„  xt>J  •  nt  '  (yS  *x  '  Vt'J' 

'  V  ’  \  ’  "V 

5y  -  -*nJ  .  ny  ■ 

J'1  ■  *5  yn  -  y5  xn  ‘ 


(1.6) 


In  general,  the  metrices  of  Eq.  (1.6)  are.  not  known  analytically  and  must  be 
determined  numerically  at  each  step  of  integration  procedure. 


The  viscous  vectors  may  be  rewritten  as  follows: 

0 


R  or  S  ■  j 


(  )  T  +  (  )  T 
x  xx  y  xy 
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(  )  T  +  (  )  T 

x  xy  wy yy 
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\  ♦  <  ♦ 


UT  +  VT 
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where  (  )  *  £  for  R  and  (  )  -  n  for  S.  ■  And 


T  ~  T .  . 

yy-  M 


■r— K£  T1  +  U  T  +  V  T 
Pr  y  'y  xy  yy 


Txx  •  T  11  (un  \  *  u5  E*>  -  f  “  <ny  v„  +  Ey  V  -  \  v  y 


t  =  v(5  ur  +  n  u_  +  £  v_  +  n  v  ) 
xy  y  £  y  n  x  £  x  n 


Tyy  ■  f  u<5y  v£  +  \  V  *  !  u<5x  "5  +  \  un  +  7  > 

■  ■  3  u5  +  \  un  +  Ey  v£  +  "y  V  +  1  v  7 

In  Eqs.  (1.7)  and  (1.2),  the  Stokes  hypothesis,  i.e.  3A  +  2y  =  0  was  used. 
Equations  (1.7)  and  (1.8)  may  be  rewritten  by  separating  terms  with  flow 
variables  differentiating  with  respect  to  £  and  ri  as  follows: 

R  =  Rx  +  R2 
S  =  Sx  +  S2 


T  -  Ti  +  T2 


where 


R  -  — 
K1  J 


b2ut  +  b3v£  -  c65kv 
b3u£  +  b4v£  -  c6’yv 


b5£c  +  (b-u  +  b-v)uf  +  (b^u  +  bAv)vr  -  ctv(u£„  +  v£y) 
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R2  or  S2  "  1 


V()  +  b9  vo 


b10u()  +  bllv() 


b7e()  +  (b8u  +  bl0V>UO  +  (b9  u  +  bllv)v() 


where  (  )  *  0  for  R2  and  (  )  =  £  for  S2. 
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where 


C2  ‘  •*  <3  \  +  ny'> 
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Equations  (1.5)  are  the  governing  equations  to  be  solved  with  the  boundary  and 
initial  conditions  described  in  section  4.  When  thin  layer  approximation  is  made 
i.e.  _9_  (  )  «  _3_  (  )  for  all  diffusion  terms,  then  one  sets  Ri  *  R2  “  S2  =  T2  “  1 
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CHAPTER  2 


NUMERICAL  ALGORITHM 


The  numerical  algorithm  described  In  this  section  Is  based  on  the  work  of 
Ref.  3.  It  Is  Intended  to  solve  the  full  Navler-Stokes  equations,  Eq.  (1.5). 
However,  the  thin- layer  approximation  will  be  made  at  the  end. 

Rewritten  Eq.  (1.5)  as  follows: 

“x  ♦  h  *  Fn  +  H  ■  r  [r1<u’ue>5  +  R2<u’Vc  +  sl<u‘un>n 

e 

♦  s2(u,uc)n  +  +  t2(u,uc)  ]  (2.: 

A  single  step  temporal  scheme  for  advancing  the  solution  of  Eq.  (2.1)  is 
...n  _  a'Ax  3  ...n  .  Ax  3  „n  .  b'  . 

^  rrr  3?au  *rrr  i?u  +rrFriU  + 


0  f(a'  •  j  *  b')  Ax2  +  Ax2] 


where  l f1  ■  U(nAt)  and  AU*1  ■  -  U°.  Substituting  UT  from  Eq.  (2.1)  Into 

Eq.  (2.2)  one  obtains 


,n  _  a 


atJ  ♦  at2  -j 
R  J 


I  n 

■  TTV  AU"'1  •>')  At2  ♦  Ax3] 


NSWC  TR  82-286 


where  En+^  ■  E(Un+^) ,  AE  =  En+^  -  E°.  A  local  linearization  can  be  achieved 
by  the  Taylor  series  expansion: 


or 


AEn  =  AnAUn  +  0  (At2). 


where  A° 


E  n 

(— )  Is  the  Jacobian  matrix.  Similarly, 


AFn  =  Bn  AUn  +  0  (At2) 


AHn  *  Kn  AUn  +  0  (At2) 


‘(w)  AU" +  Guf)AU"  +  0 1*-*2) 


-  L  n  All"  +  Pn  AU]J  +  0  (At£) 


where 


»  Ln  AUn  +  (PAU)J  +  0  (At2). 


where 


ASj  -  Hn  AUn  +  (QAU)JJ  +  0  (At2) 


ATf  «  nJ  AUn  ♦  (W^UjJ  +  0  (At2) 


where 
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AT2  -  Ng  AUn  +  (W2AU)J  +  0  (At2) 


(2.4g) 


where 


JOi.*  \\  .(OiY 

\lir  ’"2  \13U?  j 


The  cross-derivative  terms  are  treated  explicitly 

ARj  -  AR”'1  +  0  (At2) 


(2.5a) 


ASj  -  AS”"1  +  0  (At2) 


(2.5b) 


All  the  Jacobian  matrices  are  given  in  Appendix  A.  In  the  analytical 

derivation  of  Jacobian  matrices  for  viscous  portion,  it  is  assumed  that 
transport  coefficients  are  locally  constant,  i.e.  ■  ic^  ■  »  0. 

By  substituting  Eqs.  (2.4)  and  (2.5)  into  (2.3),  one  obtains  the 
spatially  factored  form: 


*■'  »_*  ■/  <  *  ■  .  •*. 
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*f?r  K>r’ + ‘“iCJtSt  “/”1 


AU" 


+  0 


[■•• 


"  ? 


1  b')  At2,  (a"-  a)  At2,  At3] 


(2.6) 


where  a"  has  been  introduced  in  the  coefficient  of  the  cross-derivative  terms 
for  notation  convenience.  3  For  second-order-accurate  schemes,  a"  should 
be  set  equal  to  a'.  However,  for  the  first-order-accurate  scheme  (a'  ^  %  +  b') 
it  is  consistant  to  set  a"  equal  to  zero.  Now,  Eq.  (2.6)  has  the  same 
temporal  accuracy  as  Eq.  (2.3)  but  is  linear  in  AU11.  In  practice,  Eq.  (2.6) 
is  implemented  by  the  sequence 


{'  [(»  - 1  )"5  -  k '  t V']  }“*■ 


RHS  of  Eq.  (2.6) 


(2.7a) 


{* ♦ ffv  [(»  - k)l  •  k q"n + ""  •  k  (Hi  ‘ vnl * au" '  “**  <2'7b> 


Prom  now  on,  it  is  assumed  that  the  thin-layer  approximation  is  applicable. 
Also  the  first-order-accurate  Euler  implicit  scheme  (a1  ■  1,  a"  “  0  and 
b'  ■  0)  is  chosen  for  the  time  integration.  As  described  in  Refs.  3  and  6, 
it  is  necessary  to  add  the  fourth-order  explicit  dissipation  terms  in  order 
to  damp  high-frequency  growth  and  thus  serve  to  control  nonlinear  instability. 
Also  the  addition  of  the  second-order  implicit  dissipation  terms  will  extend 
the  linear  stability  bound  of  the  fourth  order  terms.  Therefore,  the  final 
form  of  Eq.  (2.7)  for  thin- layer  approximation  of  the  time-dependent  compressible 
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Navier- Stokes  equations  using  Euler  implicit  time  differencing  scheme  may  be 
written  as  follows: 


I  +  AtA^  - 


<J-y5J)"]  AU*n  -  -  AX  [l»  ♦  ♦  *  -  j^.  (Sjn  ♦  T,)] 

•  **  [J'1(Vn>2  Ju]"  -  £b[J',<VS>2  Ju]"-  «•«*> 


['  t4T(B‘0"  '  k  C  +  l'"  -  k  (Nl  *  Hln>n  -  £i<j"\v>"]  4M" 


(2.8b) 


The  spatial  derivatives  appearing  in  Eqs.  (2.8)  or  (2.7)  are  approximated 
by  three-point  second-order-accurate  finite  difference 


(Vj,k*  2 AC  (fj+l »k  -  fj-l ,k) 


(2.9a) 


Wj.k  *  <fK>j.k  *  (Vl.k  ’  2fJ.k  *  fJ->.k>  (2.9b) 

(V?)2  fJ»k  *  ^4  (fj+2,k  -  4fj+l.k  *  6fj,k  '  4fJ-l,k  +  fj-2,k*  (2-9=) 

With  the  finite  difference  approximation  of  Eq.  (2.9),  a  block-tridiagonal 
system  of  the  differenced  equations  is  formed.  It  should  be  noted  that 
Eq.  (2.9c)  is  applied  to  the  RHS  of  Eq.  (2.8a)  only,  therefore  will  not 
affect  the  block  tridiagonal  system  (use  parabolic  extrapolation  for  j  *  2 
and  (J-l)).  For  the  £  -  sweep,  one  obtains 

♦j.!  4U*j_l  +  ri  K  +  Vl  4U)+1  '  (RHS)j’  3  '  (2'10) 

where  the  (BHS)^  is  a  column  matrix  which  are  computed  with  known  flow 
variables  over  the  entire  grid  points  and  and  ^  are  4x4  (m  x  1)  matrix. 

As  described  in  the  next  section,  the  boundary  conditions  at  the  axis  of 
symmetry  and  the  outflow  plane  are  imposed  implictly.  This  is  done  by 


P.  •  .  1  .  *  •  *  i  *  »  *  •  *  •  ■,  *  ■  »t  •  •  ■  ». 


v*.a 
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% 


modifying  the  coefficient  in  the  matrix  at  j  -  2  and  j  *  J  -  1,  thus 
Eq.  (2.10)  will  produce  a  bock  tridiagonal  system  as  follows: 


.  “ 

“ 

r2  *3 

au£ 

(Rhs)2 

♦2  r3  *4 
***, 

•  *  • 

•  •  • 

•  *  » 

•  •  • 

•  •  • 

AU3 
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• 

• 

• 

(Rhs)3 

•  •  • 

•  •  • 

•  •  • 

•  •  • 

•  •  • 
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*J-3  rj-2  ^J-l 

• 

• 

• 

• 

• 

<2 

(RHS)j_2 

<*J-2-*J>  rJ-i 

4UJ-1 

1 

(RHS)j.i 

•  mt 

L  J 

where,  for  Eq.  (2.8), 

At  .n  -1 

r2  -  I  +  -y-  A2  -  ex  (J2  Jj-2)  I  for  m  =  1,  2,  4 
At  n  -1 

*  *  — 2~  A 2  +  eI  ^2  for  m  *  3 

Ti  -  (1  +  2eT)  I,  j  *  3. . . .  J-2 

rj-i  "  1  +AtAj  -  2£i(Jj-i  Jj"1)n  1 

$j-l  =  “  ^Aj'l  "  GI  (Jj 

^j+i  “  T  Aj+i  "  ei  ^j+i^  1*  J*2**'*J-2 


Equation  (3.11)  applies  to  k  from  2  to  K-l.  The  solution  of  the  block 
tridiagonal  system  is  obtained  by  the  non-pivoted  LII  decomposition  method. 

Once  All*  are  obtained,  it  is  ready  to  perform  the  n-sweep.  As  described 
in  the  next  section,  the  boundary  conditions  at  the  shock  and  the  body  surface 
are  imposed  explicitly,  thus  the  flow  variables  in  these  two  boundaries 
(k  *  1  and  K)  are  updated  first  and  treated  as  known.  The  block  tridiagonal 
system  is 
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r2  *3 

r;  v. 


*K-3  rK-2  *K-1 
*K-2  rK-l 


AU? 

2 

AU? 


AU 


k-2 


AU! 


k-1 


AU. 


AU, 


AU, 


K-2 


AU 


k-i 


(2.12) 


where 


K  •  o  ♦  2ei>  1  +  >5  *  r«  -N,"k  ) .  *  -  2 ....  k-i 

♦k-i  -fK-i  -  -  h  K'  Wn «. 

K+l  =  21[Bk+l  '  5"(Mk+l*Q k+1 +  Wlk+1)]  '  EI  Jk+l)n  *»  k’2 


**  *  ~ 
AU2  =  AU2  + 


**  *  . 
“k-I  *  ^-1  -  *1 


-  «r  (o'1  v"  > 

*1  ■  T-  bk  )  ■  Ei  (JK-1  V"  1 


Equation  (2.12)  applies  to  j  *  2  to  J  -  1.  The  same  block  trldlagonal  solver 
Is  used  to  solve  (2.12).  This  completes  one  full  Integration  per  time  step. 
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CHAPTER  3 

BOUNDARY  AND  INITIAL  CONDITIONS 

As  shown  in  Fig.  1,  one  would  like  to  compute  the  flowfield  enclosed  by 
the  four  boundaries  A,B,C  and  D,  where  A  is  the  bow  shock,  B  is  the  body 
surface,  C  is  the  axis  of  symmetry  and  D  is  the  outflow  boundary.  In 
implementing  the  boundary  conditions,  one  intuitively  expects  implicit 
boundary  conditions,  to.  be  more  stable  than  explicit  one.  However,  according 
to  many  authors  13, 14, 15  this  has  not  been  their  experience.  Since  treating 
the  boundary  conditions  explicitly  is  far  more  simpler  to  implement  than  to 
do  it  implicitly,  the  shock  points  and  body  points  boundary  conditions  are 
imposed  explicitly  according  to  the  method  described  by  Kutler^  and  are 
briefly  given  in  the  following. 

3.1  Shock  Points 

The  flow  in  the  vicinity  of  the  bow  shock  is  assumed  to  be  inviscid  and 
the  Rankin e-Hugoniot  relations  are  satisfied  or  the  shock  is  tracked.  Since 
the  final  location  of  the  shock  must  come  from  the  solution,  so  it  is  allowed 
to  move.  A  quas^  steady  propagation  of  the  shock  is  assumed.  The  pressure 
behind  the  shock  is  first  determined  by  integrating  the  energy  equation  in 
nonconservative  form  as  follows: 

PT  -  -5p5  -  5pn  -  +  vtCy  +  unnx  +  vnny  +  ( 

where 

u=E+£u  +  5v» 

t  ^x  y 


v  =  ri+nu  +  r|v 
t  x  y  . 

Because  explicit  method  is  used,  the  time  step  Axs  for  shock  integration, 
Eq.  (3.1)  is  restricted  by  the  CFL  condition  (CN  =  1) ,  or 
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where  for  A:  kQ=  £t,  ^  -  £x ,  k2  »  £y;  and  for  B:  kQ  =  nt»  \  =  Ux»  k2  ■  Uy; 

and  the  constant  0.9  is  a  safety  factor  which  must  be  less  than  1.  It  should 
be  noted  that  Axs  is  different  from  the  At  used  in  the  integration  of  the 
interior  points,  this  means  that  the  calculation  cannot  be  time-accurate. 
However,  this  does  not  prevent  one  to  obtain  steady  state  solution. 

Knowing  the  pressure,  the  shock  velocity  can  be  determined  as  follows 
(see  Fig.  2) : 


FIGURE  2.  NOTATION  FOR  SHOCK  POINT  BOUNDARY  CONDITION 


where 


a  M  -  u. 
»  x  1 


(3.4) 


M. 


i*[2 


(Y  +  1)  +  (Y  -  1) 


])! 


=  q^cosS 

a®  -  orPoo/pJ*5 

The  shock  angle  6  is  a  function  of  the  metrics 


6 


tan_1(-Ti  /r|  )  .  . 

y  x  shock 


(3.5) 
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where  r)y  and  r|x  are  defined  in  Eq.  (1.6).  From  the  Rankine-Hugoniot  relations 
the  density  behind  the  shock  can  be  determined  by 


/(-£•£) 


The  velocity  components  behind  the  shock  in  cylindrical  coordinates  are: 

U2  =  q^sin  6  +  locos'? 


(3.6) 


(3.7) 


=  q^sinficosS  - 


where 


u2  -  2(1  -  Mx“)  am  /  (Y  +  1)M, 


[  1  +  un . 
xj  1 


The  value  of  e  can  be  obtained  from  Eq.  (2.3a)  with  known  P2,  P2>  V2* 

Once  the  shock  velocity  qs  of  Eq.  (4.4)  is  known,  the  new  shock  position  at 
time  T  +  At  can  be  determined  by  propagating  the  shock  along  a  (  =  constant 
line,  or  equivalently  in  the  n  direction  with  a  velocity. 


^s^  =  qs/cos8’ 


(3.8) 


where  8'  =  0  -  6  and  0  =  tan  ^  (£x/£y).  This  determines  the  new  x  and  y  values 
of  the  shock  points  and  subsequently  the  new  x  and  y  values  of  the  interior  points. 

3.2  Body  Points 

The  tangency  condition  for  inviscid  flow  at  the  body  surface  is  imposed 
explicitly  using  Kentzer's^  scheme  through  characteristic  analysis  (see  Appendix 
B)  to  yield  the  following  equation  for  the  pressure  time  derivative  on  the  body: 

PT  =  "  /nP2+n  2  1<VT  +  u(Vt  +  v(Vt  ]  +  3  ^x2+ny2  pn  _  SPS 


x  y 


-pa2  (f  ur  +  ?  v,  +  ri  u  +  nv  +  -  ) 
■X.  I  y  K  xn  'y  n  y  ’ 


where 


+  g 

2+n  2 
x  y 


[nx(pu  u^  +  5xp^)  +  ny(puv^  +  Syp^)] 


(3.9) 


S  -  5t  +  +  v5, 


v  =  nt  +  unx  +  vny  =  o 


^Kentzer,  C.  P.,  "Discretization  of  Boundary  Conditions  on  Moving  Discontinuities," 
Lecture  Notes  in  Physics,  published  by  Springer-Verlag,  No.  8,  Sep  1970,  pp.  108-113. 

26 


NSWC  TR  82-286 


In  the  numerical  code,  the  time  derivative  terms  were  neglected  and  the  time  step 
was  modified  by  a  safety  factor  of  0.2  for  stability  reason.  The  stagnation 
pressure  was  enforced  at  the  stagnation  region  (J  *  1  and  2)  as  follows: 


A"  (P4  rfP3  l+f  Pt)/6.25 
B  '  *<»3,1  -  *t  -  T  A)/’ 


^  (3.10) 


J 


where  R  is  the  prcssu'e  at  point  (J,K)  and  Pt  is  the  total  stagnation  pressure. 
It  was  r6und,  Eq.  3.10  nolps  to  speed  up  the  convergence  process. 

Once  the  surface  pressure  is  determined,  the  remaining  variables  on  the  body 
surface  are  determined  by  the  following  isentropic  relations: 


where 


Pb  -  (Pb/s)1/Y 


Ub  "  %  COS  0b 
Vb  =  qb  Sin  9b 


>  (3.11) 


»  h.  j[2YM»2  -  (Y-1)]/(Y+1)  Y  _P® 
’  P1Y  i(Y+l)H»2/[(Y-l)M»2+2]  ’  PJ 


is  the  stagnation  entropy,  Pt,  pt  are  the  stagnation  pressure  and  density. 


[%(Y+1)M»21  y~1 


6  {[2YH*2  -  (Y-1)]/Y+1)1^(Y"1)  °° 


(3.13) 


rt  'fe) 


and  0^  is  given  by 


w""1(-Wb0dy 


(3.14) 
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with  nx  and  defined  by  Eq.  1.6. 

To  simulate  viscous  flows,  the  no-slip  boundary  condition  requires 

u  ■  v  *  0 

when  the  body  surface  is  not  moving,  Eq.  (3.15)  also  implies  u  -  v  -  0.  To 
determine  the  surface  pressure,  it  is  assumed  that  the  normal  pressure 
gradient  over  the  first  3  grid  points  above  the  body  surface  is  zero,  or 


p„  “  — — 5T  r<5  n  +  K  n.  )Pc  +  (n  2  +  n  2)p„1  =  o  (3.1 
n  (n  2  +  n  2)*5  L  x  x  yy  ?  x  y  nj 

x  y 

Similarly,  for  an  adiabatic  wall  boundary  condition,  the  temperature  may  be 
determined  by 

(£  n  +  ?  n  )t_  +  (n  2  +  n  2)t  =  o  (3.1 

x  x  *y  y  5  x  y  ri  v 

The  £  and  n  derivatives  in  Eqs.  (3.16)  or  (3.17)  are  differenced  using  a 
second-order  central  difference  formula  for  the  ^-derivatives  and  a  three- 
point  one-sided  formula  for  the  n-derivatives .  This  results  in  a  tri¬ 
diagonal  system  of  equations  which  can  be  solved  to  yield  the  pressure  and 
temperature.  For  a  constant  temperature  wall,  the  temperature  along  the  wall  is 
kept  constant  at  its  initialized  value  throughout  the  entire  convergence 
process.  Once  pressure  and  temperature  are  known,  the  density  is  determined 
from  the  equation  of  state. 

3.3  Plane  of  Symmetry 

The  axis  of  symmetry  line  is  bypassed  by  choosing  the  first  two  £  = 
constant  lines  to  straddle  the  axis  or  stagnation  streamline.  The  plane  of 
symmetry  boundary  is  then  enforced  by  the  reflection  principle.  The  flow 
variables  are  either  even  or  odd  functions  with  respect  to  the  plane  of 
symmetry,  or 


P(l»k)  -  p(2,k),  v(l,k)  =  -v(2,k) 

u(l,k)  *  u(2,k) ,  e(l,k)  =  e(2,k) 

The  boundary  condition  at  the  plane  of  symmetry  is  imposed  implicitly. 
3.4  Outflow  Points 


(3.1: 


Since  velocity  at  majority  of  the  grids  in  this  plane  is  supersonic,  a 
simple  linear  extrapolation  of  the  conservative  variables  is  used.  For  those 
points  near  the  surface,  the  flow  is  subsonic  there,  thus  some  error  is  intro¬ 
duced.  The  supersonic  outflow  boundary  condition  is  imposed  implicitly,  or 


Q(J,k)  -  2Q(J-l,k)  -  Q(J-2,k) 


(3.1 
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3.5  Initial  Conditions 

To  start  the  calculation,  an  initial  flowfield  must  be  provided.  In  the 
numerical  code,  two  starting  methods  are  provided:  (1)  starting  the  calculation 
with  a  sphere  at  a  given  free  stream  Mach  number  or  (2)  reading  in,  point  by 
point,  the  shock  and  body  locations  along  rays  (£  *  constant  lines).  For  the 
indented  nosetip  shapes  described  in  this  report,  the  calculation  always  starts 
from  a  sphere  and  let  the  body  gradually  deformed  to  the  desired  shape  by  giving 
a  set  of  values  for  6  along  each  ray  (see  Fig.  1).  For  sphere,  a  good  guess  for 
the  shape  and  position  of  the  shock  is  made  based  on  known  spherical  blunt  body 
solutions  as  follows: 


SQ  -  0.78  { [ (Y-1)M»2  +  2]/[(Y-l)M»2]> 
yb  -  2.376  -  0.1834Moo  +  0.01036M»2,  3  <  <  10 

-  1.576  -  0. 001804*  -  10),  >  10 

ps  *  yt,/t2(1+5o>l 

y8  -  2P8(x8  +  60)  (3.20 

where  6q  is  the  stand  off  distance  at  axis  of  symmetry,  y.  and  P8  are  empirical 
constants  as  a  function  of  Mach  number  and  yg  and  x8  are  the  shock  position.  As 
one  can  see  a  parabolic  shock  shape  is  assumed. 

Once  the  shock  and  body  locations  are  determined,  the  flow  variables  along 
the  shock  are  obtained  by  assuming  a  zero  shock  velocity  and  applying  the 
Ranklne-Hugonlot  relations.  On  the  body,  a  Newtonian  pressure  distribution  and 
isentropic  relations  are  applied  to  provide  the  flow  variables.  The  equations 
to  implement  the  initial  flow  variable  on  shock  and  body  are  listed  in  the 
following: 

On  shock: 


P8  -  (  [2Y(Mb>  sin  0S)2  -  (Y-l)  ]/ (Y+l)  }  P*, 

P8  -  [(Y-l)Ofcosin  68)2/[  (Y-l)  (Moo  sin  0)2  +  2] }  p*, 

u8  -  (1-2(M»2  sin2  08  -  1)/[(Y+1)M»2])  q*, 

v8  -  (2 (Moo2  sin2  08-l)cos  0S/[(Y+1)M»2  sin  9S]  )qco  . 


(3.21 
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At  stagnation  region 


s  -  H2^2-  P’-Dim+i)  j Y  Poo 

/  (Y+D^/fCY-DH*,^]  I  Jv 

00 


t  { [2YH»2-  (Y-l)  ]  /  (Y+l)  >1//y_1  ftCY+DM-  ]  Poo 

t/Pt  =  [1 +%(Y-1)^2]  • 

Kqo 

On  body: 

Pv  -  Pootfcr-  -  1)(1.0  -  1.02  sin20  +  0.12  sin40)  +  1.] 

D  Poo 


vb  ’  %  8in0b 


pb  - 

(Pb/S)*8 

I 

qb- 

\  2Y  ,  pt 

Pb 

/*\ 

1 

I—* 

T>  1 
rt  | 

1 

Pb 

ub " 

1  %  CO80b  | 

The  flow  variables  at  each  nodal  point  within  the  shock  layer  are  linearly 
interpolated.  In  order  to  maintain  a  constant  total  enthalpy  for  the  initial 
flowfield  at  each  nodal  point,  a  modification  of  the  interpolated  valocity 
components  is  made 


where 


ui  “  ui 
vi  *  vi  ^q^q^ 


(3 


q  -  [^-/Pt/Pt  -  ^/] 

q  ■  +  v'2 
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CHAPTER  4 


INDENTED  NOSETIP  SHAPES  AND  COMPUTATIONAL  MESH 


4.1  Indented  Nosetip  Shapes 

In  this  report,  the  series  of  Indented  nosetips  as  shown  In  Figure  3  are  of 
primarily  concerned.  These  shapes  are  composed  of  arcs  and  straight  lines.  Thr 
portions  can  be  identified:  (I)  expansion  corner  near  the  flat  nose,  (II) 
compression  turn  and  (III)  expansion  shoulder.  The  detail  description  of  the 
models  is  given  in  Table  1.  It  should  be  emphasized  that  these  nosetips  present 
serious  computational  difficulties  because  of  the  sharp  comers  and  convex  and 
concave  curvatures.  Experimental  measurements  of  bow  shock  location  and  surface 
pressure  for  model  1  to  3  are  reported  in  Ref.  7.  Also  velocity  and  density 
data  have  been  obtained  in  Ref.  8  for  model  4.  Numerical  calculations  will  be 
compared  to  these  results. 


Table  1  Description  of  Nosetip  Models 


Mod«l  Para¬ 
meter 


I  UIJUUU3JUUUUU 

X/L 


Subscript  n 


XPn/L  0.053  0.077  0.579  1.508  2.474  3.140  5.906 
YP-/L  0.676  0.715  1.101  2.037  3.274  3.651  3.990 
ftn/L  4.320  0.062  5.389  1.000 
52.00  38.00  7.000 


XPn/L  0.038  0.069  0.572  1.475  2.712  3.378  6.144 

YP./l  0.569  0.615  0.906  1.692  3.274  3.651  3.990 

Rn/L  4.320  0.062  3.137  1.000 

fate?  60.00  38.00  7.000 


0.025  0.064  0.564  1.454  2.929  3.595  7.631 

0.463  0.514  0.717  1.388  3.274  3.651  3.990 

4.320  0.062  2.154  1.000 

68.00  38.00  7.000  * 


XPn/L  0.  0.263  0.830  1.693  2.099  2.506  5.225 

YPn/L  0.266  0.642  0.848  1.605  2.386  2.S69  3.000 
K./L  •  0.400  1.585  0.533 

70.00  27.50  7.000 


FIGURE  3.  NOSETIP  SHAPES 


4.2  Computational  Grids 

As  shown  in  Figure  1,  the  computational  grids  are  formed  with  constant 
spacing  rays.  The  first  two  rays  (j«l  and  2)  must  straddle  the  axis  of 
symmetry.  Hence 

A8  «  8  / (J-l. 5) 

max 


8  -  (j-1. 5)*A8 
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To  deform  a  sphere  to  the  indented  nosetip  shapes,  the  distance  of  deformation  6(6) 
must  be  calculated  for  each  ray.  Expressions  for  6(6),  x(6) ,  y(8)  are  given  in 
Appendix  C. 

The  distribution  of  points  along  a  ray  between  the  shock  (x8,y8)  and  the 
body  (xb,yb)  is  given  by 

I(k)  -  1  +  B[l-<f^)1"b]/[H(|i|)1“b]  (4.2a) 

where 


b  -  (k-1) / (K-l)  (4.2b) 

and  6  is  a  free  parameter.  A  uniform  distribution  of  grids  is  obtained  as  8  -+•  » 
and  a  strongly  clustering  of  points  near  the  body  surface  is  obtained  as  8  -►  1. 
Table  2  gives  the  first  6  values  of  5(k) ,  k-l, 7  for  8  -  1.01,  1.005  and  1.001. 

TABLE  2.  MESHES  USED  FOR  HEMISPHERE- CYLINDER  AND  SPHERE-CONE  CALCULATION 


Mesh 

8 

Points 

O.lYs 

First  6  points  | 

nj/Ys 

n2/Ys 

nj/Ys 

n4/Ys 

ns/Ys 

n6/Ys 

A 

1.01 

12 

.0025 

.0055 

.0093 

.014 

.0199 

B 

1.005 

17 

.0011 

.0024 

.0039 

.0058 

.0081 

n 

1.001 

20 

.0003 

.0011 

.0017 

.0024 

■ 

1.005 

17 

.0011* 

.0024 

.0039 

.0058 

.0081 

The  distribution  of  a(k)  once  initiated  is  kept  unchanged  throughout  the 

calculation  even  though  the  grid  points  may  vary  at  each  time  step  as  the  shock 

or  body  location  varies.  For  invlscid  calculation,  uniform  distribution  of 
point  across  the  shock  layer  is  always  used  with  8  -  10$  in  Eq.  (4.2).  The 
clustering  of  points  near  the  surface  must  be  used  for  viscous  calculations. 

4.3  Criterion  for  Convergence.  The  convergence  of  solution  is  judged  by  the 
convergence  of  shock  speed.  From  experience  for  smooth  nosetip,  when  the 
non-dimensional  shock  speed  reaches  10“ 3  to  10”5 ,  all  the  flow  variables  remain 
essentially  constant  and  the  residue  of  the  solution,  which  is  given  by  the  RMS  of 
Eq.  (2.7),  is  in  the  order  of  10~5  to  lO”?  over  all  the  grid  points.  The  shock 
speed  does  not  converge  monotonously,  it  will  oscillate  but  decrease  in  the  averaged 
magnitude.  For  all  cases  computed  for  smooth  nosetlps,  a  converged  solution  may 
be  obtained  in  about  300  to  1000  time  steps  depending  on  the  mesh  selected. 

In  general,  the  convergence  rate  slows  down  rapidly  for  a  strongly  clustered  mesh 

such  as  Mesh  C;  see  Table  3  for  the  hemisphere-cylinder  case. 
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TABLE  3.  CONVERGENCE  RATE  FOR  HEMISPHERE-CYLINDER 


Mesh 

Number  of  Time 
Integrations 

Shock  Speed 

Remark 

A 

300 

.51  x  10"4 

CM  - 
eE  " 
eI  ~ 

75  and  final 

0.02  and 

3  Eg  for  all  cases. 

B 

600 

.68  x  10"4 

Mesh 

eg  ■ 

C  starts  with 

0.1  for  400  steps. 

C 

1000 

.16  x  10"2 

For  Indented  nosetips,  the  convergence  of  shock  speed  becomes  very  slow  when 
It  reaches  the  order  of  0.01  with  long  oscillation  cycle  in  both  inviscid  and 
viscous  calculations.  Under  such  conditions,  the  flow  variables  near  the  body 
are  essentially  unchanged  for  over  a  few  hundred  time  steps  and  calculations  are 
then  considered  completed.  The  residue  may  go  up  to  10“  ^  in  some  grid  points. 

4.4  Time  steps 


The  time  step  At  used  for  each  time  integration  is  determined  from  the  input 
Courant  number  CN  according  to 


At 


CN 

amax 


(4.5) 


where  amax  Is  given  by  Eq.  (3.3)  and  is  the  maximum  of  the  eigenvalues  of  the 
matrices  A  and  B  over  all  the  interior  points.  For  inviscid  flow  calculations, 
the  allowable  value  of  CN  is  2  for  both  smooth  and  indented  nosetip  calculation. 
For  viscous  flow  calculation,  the  allowable  Courant  number  is  problem  dependent 
and  also  related  to  the  dissipation  coefficients  used.  For  smooth  nosetip  laminar 
flow  calculations  eg.  ■  0.02  and  ej  =  3£g  are  used  and  a  CN  of  75  has  been  achieved 
(can  also  run  for  CN  *  150  if  let  Eg  ■  0.1).  In  some  situations,  a  larger  value 
of  eg  is  required  at  the  beginning  of  the  convergence  process  and  may  be  reduced 
afterward  as  the  flowfield  gradually  approaching  the  final  solution.  For 
indented  nosetip  laminar  flow  calculation,  CN  has  been  reduced  to  50  or  25  in  some 
part  of  the  calculation.  The  CPU  time  required  for  the  integration  is  about 
0.000851  per  time  step  per  mesh  point  using  GDC-7600  computer. 


4.5  Dissipation  Coefficients 

Experiments  about  the  effect  of  adding  dissipation  terms  on  the  solution  were 
conducted  for  the  cases  of  sphere  and  hemisphere-cylinder  for  both  inviscid  and 
laminar  flow  calculations.  For  these  smooth  noses,  no  significant  difference 
between  the  solutions  with  Eg  =  0.02  and  0.001  (£j  =  3sg)  can  be  found.  However, 
if  let  Eg  *  Ej  ■  0  the  solution  diverges.  Thus,  for  smooth  nosetip  Eg  ■  0.02  is 
used  for  all  calculations.  For  indented  nosetip  calculations,  it  is  necessary 
to  have  a  larger  value  of  Eg  in  the  order  of  0.2  to  0.4  in  order  to  have  a 
solution.  No  assessment  about  the  influence  of  the  dissipation  terms  on  the 
solution  can  be  made  at  present  time.  How  good  the  solution  is  can  only  be  judged 
by  comparison  to  experimental  data  for  the  indented  nosetips.  More  discussion 
about  the  effect  of  dissipation  are  provided  at  the  presentation  of  results. 
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CHAPTER  5 

INVISCID  FLOW  CALCULATIONS 


5.1  Code  Verification 

A  copy  of  the  research  computer  code  described  in  Ref.  1  was  obtained  from 
NASA  Ames  Research  Center,  referred  to  as  K-C-L  code.  The  inviscid  portion 
(mainly  contributed  by  P.  Kutler)  of  the  code  is  first  verified  for  the  case  of  a 
sphere  at  =  5.96.  Grid  points  of  6  (normal  direction)  x  10  (streamwise 
direction),  12  x  19  and  12  x  32  have  been  used  and  results  obtained  for  surface 
pressure  agree  well.  The  total  number  of  time  steps  used  in  each  calculation 
is  600  and  the  final  shock  speed  (normalized  by  free  stream  sound  speed)  reached 
is  ih  the  order  of  10”^  to  10“®  which  is  considered  to  be  sufficiently  small  for 
the  results  to  be  judged  as  steady  state.  The  maximum  total  enthalpy  error  in  the 
flowfield  varies  from  0.7%  for  the  6  x  10  grid  to  0.1%  for  the  12  x  32  grid;  this 
relation  is  almost  linearly  proportional  to  the  grid  points  used.  Therefore,  the 
code  seems  to  behave  consistently  for  sphere.  Next,  the  effects  of  explicit 
dissipation  terms  used  in  the  numerical  scheme  is  examined.  The  value  of  eE  is 
varied  from  0.4  to  0.005  and  the  resulting  surface  pressure  agree  to  the  first 
two  digits.  When  eE  ■  0  or  ^  0.6,  the  calculation  diverges.  Thus,  the  effect 
of  explicit  dissipation  on  the  calculated  results  for  flow  over  sphere  is  small. 

As  will  be  discussed  later,  they  are  not  so  small  for  the  indented  nosetip.  In 
the  above  calculations,  the  maximum  Courant  number  is  2.5.  The  results  obtained 
with  12  x  32  grid  are  compared  to  the  work  of  Inouye  and  Lomaxl?  and  the  measured 
data  of  Baerl®  for  the  surface  pressure  as  shown  in  Fig.  4a  and  to  the  measured 
data  of  Sedney  and  Kahili  for  the  density  distribution  as  shown  in  Fig.  4b.  The 
agreements  are  seen  to  be  satisfactory  except  for  p/Poo  *  5.2  in  Fig.  4b.  It 
should  be  pointed  out  that  the  calculated  result  for  p/poo  «*  5.2  agrees  well  with 
the  numerical  results  of  Shubin  et  al*  who  solve  the  steady  Euler  equations  in 
conservation  law  form  with  an  entirely  different  numerical  method. 


^Inonye,  M. ,  and  Lomax,  H.f  "Comparison  of  Experimental  and  Numerical  Results  for 
the  Flow  of  a  Perfect  Gas  About  Blunt-Nosed  Bodies,"  NASA  TN  D-1426,  Sep  1962. 
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Baer,  A.  L.,  "Pressure  Distribution  on  a  Hemphisphere  Cylinder  at  Supersonic  and 
Hypersonic  Mach  Numbers,"  AEDC  TN  61-96,  Arnold  Engineering  Development  Center, 
1961. 


^Sedney,  R. ,  and  Kahl,  G.  D.,  "Interferometric  Study  of  the  Blunt  Body  Problem," 
Ballistic  Research  Laboratory  Report  No.  1100,  1960. 
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a.  SURFACE  PRESSURE 


FIGURE  4.  COMPARISON  OF  SURFACE  PRESSURE  AND  DENSITY  FIELD  BETWEEN  CALCULATION 
AND  EXPERIMENT 

5.2  Special  Calculation  Procedure  for  Indented  Nosetips 

Difficulties  arise  when  apply  the  code  to  compute  the  flowfield  around  the 
indented  nosetips  given  in  Fig.  3.  In  order  to  represent  the  nosetip  shape 
reasonably  well,  a  sufficient  number  of  grid  points  must  be  located  in  the  areas 
where  the  body  geometry  changes  rapidly.  The  code  fails  to  carry  out  the 
computation  with  both  the  starting  methods  provided.  These  two  starting  methods 
are  either  to  deform  from  a  sphere  or  to  specify  the  locations  for  the  initial  body 
and  shock.  Only  with  a  coarse  grid,  for  example  12  x  19,  can  a  solution  be 
obtained  with  the  body  not  well  represented,  particularly  near  the  nose  tip.  The 
solution  obtained  then  suffers  an  error  of  the  total  enthalpy  in  the  flowfield 
as  high  as  20  percent  for  a  severely  indented  nosetip.  Such  solution  is  obviously 
unsatisfactory. 

Experience  Indicate  that  the  problem  always  occur  at  the  shock  boundary  and 
results  in  a  shock  speed  which  continuously  increases  as  the  time  progresses  and 
finally  leads  to  the  breakdown  of  the  computation.  One  of  the  reasons  that  the 
code  behaves  in  this  manner  is  perhaps  that  the  shock  boundary  points  are 
treated  explicitly  and  in  a  quasi-steady  manner  (it  is  not  time  accurate);  thus  the 
starting  flowfield  must  be  sufficiently  close  to  the  final  flowfield  in  order  for 
the  computation  to  carry  through  smoothly.  A  fundamental  improvement  of  the 
algorithm  is  to  Incorporate  an  unsteady  shock  boundary  condition  and  to  make  the 
entire  calculation  time  accurate  with  implicit  treatment  of  shock  and  body 


boundary  conditions.  This  requires  considerable  modification  of  the  code  (in 
fact,  a  new  code)  and  is  not  persuaded  in  this  report.  The  alternative  approach 
would  be  to  proceed  slowly  by  gradually  approximating  the  body  geometry  and  this 
has  been  successfully  accomplished  as  described  in  the  following  paragraphs. 

The  optimal  procedure  for  computing  the  flow  about  an  indented  nosetip  is  to 
obtain  a  solution  with  a  coarse  mesh.  Additional  grid  points  are  then  added  in 
areas  of  rapid  geometry  variation  and  a  new  solution  is  obtained  using  the  previous 
one  as  a  starting  guess. 

The  calculation  made  for  Model  2  at  =  6  is  used  as  an  example  to  demonstrate 
the  procedure.  To  begin  with,  a  coarse  mesh  (typically  12  x  19)  for  sphere  with 
equal  angular  increment  rays  (i.e.  £  =  constant  lines,  the  first  two  rays  have  to 
straddle  the  axis)  is  used  and  a  calculation  is  performed  with  the  body  being 
guadually  deforming  in  time  along  each  ray  to  the  desired  nosetip  shape.  The 
surface  pressure  obtained  is  plotted  versus  the  arch  length  as  shown  in  Fig.  5  by 
the  broken  line.  The  intersection  point  of  the  ray  and  the  body  surface  is 
numbered  along  the  body  as  shown  in  the  sketch  in  Fig.  5.  This  solution  is  poor 
because  the  expansion  corner  is  apparently  missed  and  the  maximum  error  in  total 
enthalpy  is  11.87%.  An  obvious  cause  is  an  insufficient  number  of  grid  points  in 
the  stagnation  region  and  around  the  expansion  corner.  New  rays  are  added  which 
contain  the  same  number  of  uniformly  distributed  grid  points  between  the  body  and 
and  the  shock.  All  of  the  original  grid  points  are  kept  where  they  are.  The  flow 
variables  at  the  new  grid  point  are  interpolated  from  the  flowfield  solution  just 
obtained.  Thus  a  new  initial  condition  is  obtained  to  continue  the  calculation. 

uemMSFMDram 


FIGURE  5.  PROGRESS  OF  SURFACE  PRESSURE  AS  A  FUNCTION  OF  ADDING  GRID  POINTS 

As  shown  in  Fig.  5,  significant  change  in  surface  pressure  distribution  is 
found  after  adding  three  more  rays  (denoted  by  if 20,  21  and  22).  The  error  in  total 
enthalpy  also  reduces  to  5.09%.  The  same  procedure  is  repeated  by  adding  ray  if 23 
(the  end  point  of  the  circular  arc  of  the  expansion  corner)  and  if 24.  Now,  the 
surface  pressure  looks  like  flow  over  a  corner.  In  order  to  represent  the  comer 
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better,  the  #25  ray  is  added  in  the  middle  of  the  circular  arc.  This  addition  is 
seen  to  catch  the  minimum  surface  pressure  of  the  expansion  corner.  At  this 
point,  it  is  felt  that  the  body  geometry  is  well  represented.  A  further  attempt 
to  add  rays  between  #24  and  #7  fails  to  obtain  a  solution.  This  part  of  the 
calculation  is  most  difficult  and  time  consumming.  Each  addition  of  grid  points 
requires  about  600-1200  time  steps  to  guarantee  that  the  solution  will  converge. 

At  first,  it  seems  that  the  solution  obtained  with  12  x  35  mesh  is  sufficiently 
accurate  as  shown  in  Fig.  5.  Surprisingly  enough,  when  three  more  rays  (#26,  27 
and  28)  are  added  at  the  expansion  shoulder,  the  surface  pressure  there  changes. 
Further  addition  of  rays  in  the  downstream  of  #28  does  not  have  a  large  effect  on 
surface  pressure.  Yet  the  addition  of  rays  in  front  of  #26,  i.e.  # 29-33,  are  seen 
to  catch  the  details  of  a  process  of  recompression  and  expansion  and  recompression 
again.  A  further  increasing  the  number  of  rays  between  #29  and  #28  will  change 
the  surface  pressure  locally,  but  the  general  trend  of  the  curve  remain  the  same. 

The  addition  of  rays  between  the  compression  turn  and  the  expandion  shoulder  will 
not  present  any  difficulty  in  obtaining  a  converged  solution. 

Next,  the  effects  of  explicit  dissipation  used  in  the  code  is  examined.  This 
is  done  by  reducing  the  value  of  Eg.  As  shown  in  Fig.  5,  the  minimum  value  of 
eE  required  to  obtain  a  converged  solution  is  0.2  (or  eE/At  25)  and  the  difference 
in  surface  pressure  between  eE  =  0.2  and  0.3  is  small.  Quantitative  comparison 
of  results  with  and  without  the  explicit  dissipation  for  the  indented  nosetips 
investigated  in  this  paper  is  not  possible  because  of  the  lack  of  a  solution  with 
eE  =  0.  However,  Wardlaw*  has  calculated  a  midly  indented  nosetip  at  Moo  =  5 
using  MacCormicks ' s  explicit  method  without  explicit  dissipation.  The  K-C-L  code 
has  been  run  on  the  same  configuration  with  eE  =  0.2  and  0.005  (eE/At  1).  A 
comparison  of  calculated  surface  pressure  indicate  that  the  eE  =  0.005  run  agrees 
better  with  Wardlaw's  surface  pressure  calculation.  There  is» a  maximum  of  8% 
discrepancy  between  the  eE  =  0.2  and  0.005  solution  in  the  location  with  large 
flow  gradient. 

The  procedure  used  to  carry  out  indented  nosetip  calculation  can  be  summarized 
as  follows:  (1)  Use  a  coarse  mesh  for  a  sphere  to  start  the  calculation  and 
deform  the  sph*  re  along  each  ray  to  the  desired  nosetip  shape.  (2)  Add  new  rays 
(a  few  rays  at  a  time)  to  the  critical  areas  featuring  rapid  variation  in 
geometry  and  obtain  a  new  converged  solution.  (3)  Reduce  the  value  of  the 
explicit  dissipation  coefficient  to  the  minimum  value  producing  a  converged 
solution.  This  entire  calculation  sequence  requires  a  net  CPU  time  of  about 
15-25  minutes  on  a  CDC  7600  computer. 

5.3  Flowfield  for  Indented  Nosetips 

Following  the  calculation  steps  described  in  the  previous  section,  the 
inviscid  solutions  for  the  series  of  nosetips  given  in  Fig.  3  were  obtained  and 
results  are  present  in  the  following  paragraphs. 

5.3.1  Comparison  of  Surface  Pressure  and  Shock  Location 

The  inviscid  solution,  flowfield  pictures  and  the  measured  data  for  Model  1  to 
4  are  presented  in  Figs.  6  to  9  respectively.  The  degree  of  indentation 
increases  with  the  model  number.  Model  1  is  the  least  indented  nosetip.  As  seen 

ff 

Wardlaw,  private  communication. 


37 


NSWC  TR  82-286 


from  the  flowf ield  picture  (not  shown) ,  the  flow  remains  attached  over  the 
expansion  corner  (0.679  <  S/L  <  0.725)  and  separation  occurs  near  the  comprejsion 
turn  (1.36  <  S/L  <  2.5).  Good  agreement  of  surface  pressure  (Fig.  6b)  between 
numerical  solution  and  measured  data  up  to  S/L  ^  1.3  is  obtained.  After  the  flow 
separation,  a  sudden  increase  in  the  measured  surface  pressure  can  be  easily 
understood  by  the  concept  of  displacement  thickness  which  suddenly  thicken  the 
body  at  the  point  of  separation.  Since  the  separation  bubble  is  small,  the 
difference  between  the  calculated  and  measured  surface  pressure  is  also  small. 

The  compression  turn  is  covered  by  the  separation  bubble,  hence  the  slight  expansion 
and  recompression  before  the  expansion  shoulder  as  predicted  by  the  inviscid 
solution  is  not  given  by  the  measured  data.  It  is  difficult  to  identify  from  the 
flowf ield  picture  where  the  flow  reattaches,  but  it  seems  likely  that  the  flow 
reattaches  before  the  expansion  shoulder  (4.3  <  S/L  <  5.0)  because  of  the  good 
agreement  obtained  between  calculated  and  measured  surface  pressure  after  S/L  =  4.0. 
The  comparison  of  experimental  and  computed  shock  locations  is  in  good  agreement 
as  shown  in  Fig.  6a.  The  inviscid  solution  gives  a  thinner  shock  layer  as  expected. 
Also  shown  in  Fig.  6a  are  the  subsonic  region  and  the  sonic  lines. 

The  results  for  Model  2  are  given  in  Fig.  7.  Experimental  data  indicate  that 
the  location  of  separation  and  the  surface  pressure  are  sensitive  to  the  Reynolds 
number,  which  accounts  for  the  wide  spread  in  experimental  measurements.  The 
flow  separates  from  the  downstream  of  the  expansion  corner  (0.566  <  S/L  <  0.628). 

The  compression  turn  (1.2  <  S/L  <  2.5)  is  submerged  in  the  separation  bubble.  The 
flow  probably  reattaches  at  S/L  ^4.0  and  the  expansion  shoulder  is  located  at 
4.4  <  S/L  <  5.2.  Because  of  a  larger  separation  bubble  in  this  case,  the  difference 
between  calculated  and  measured  shock  location  (Fig.  7a)  and  surface  pressure 
(Fig.  7b)  in  the  separation  region  is  increased  also.  The  subsonic  region 
associated  with  the  compression  turn  also  increases  in  this  case. 


x/l  (a) 


FIGURE  6.  COMPARISON  OF  INVISCID  SOLUTION  WITH  EXPERIMENTS  FOR  MODEL  1  AT  JL 
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The  most  severely  indented  nosetip  of  this  group  is  Model  4.  As  shown  by  the 
holograph  of  Fig.  9a,  the  flow  is  separated  immediately  at  the  beginning  of  the 
expansion  corner  (0.266  <  S/L  <  0.75)  and  the  separation  region  extends  all  the 
way  to  the  expansion  shoulder  (3.42  <  S/L  <3.95).  The  compression  turn  is 
located  at  1.35  <  S/L  <  2.52,  and  is  fully  submerged.  Therefore,  good  agreement 
between  inviscid  solution  and  measured  data  for  surface  pressure  (Fig.  9c)  can 
only  be  found  in  the  portion  of  stagnation  region  and  for  S/L  £  3*5-  The 
comparison  of  shock  location  (Fig.  9b)  is  also  poor.  A  further  analysis  of  this 
case  using  the  concept  of  effective  body  will  be  described  later. 


FIGURE  9.  COMPARISON  OF  INVISCID  SOLUTION  WITH  EXPERIMENTS  FOR  MODEL  4  AT 
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5.3.2  Contour  of  Pressure  and  Density  and  Velocity  Vector 

Plots  of  the  inviscid  solution  for  the  pressure,  velocity  vector  and  density 
distribution  in  the  shock  layer  for  Model  1-4  are  shown  in  Fig.  10.  The  range 
and  number  of  constant  pressure  and  density  contour  are  given  in  Table  4..  From 
the  pressure  and  density  contour  plots,  it  is  seen  that  there  are  no  strong  embedded 
shock  or  slip  surfaces  in  the  shock  layer  for  all  the  models.  On  the  velocity 
plots,  it  is  interesting  to  see  that  in  the  compression  turn  area  the  velocity 
magnitude  distribution  toward  the  body  surface  first  decreases  and  then  increases 
to  form  a  retarded  velocity  regime.  As  the  degree  of  indentation  increases,  this 
retarded  velocity  region  also  increases.  For  Model  4,  part  of  the  flow  in  the  core 
of  the  retarded  velocity  region  reverse  its  direction  to  form  two  vortices.  One 
would  generally  not  expect  inviscid  flow  solution  to  behave  like  this,  but  such 
solution  should  not  be  excluded  either  since  all  the  boundary  conditions  are 
satisfied.  Whether  the  explicit  dissipation  (like  effective  viscosity)  is  solely 
responsible  for  the  resulting  velocity  field  is  not  clear.  Nevertheless,  the 
appearance  of  the  retarded  velocity  region  does  suggest  that  a  flow  separation  is 
likely  to  occur. 

TABLE  4.  VALUES  FOR  PRESSURE  AND  DENSITY  CONTOUR 
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FIGURE  10.  INVISCID  FLOWFIELD 
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(c)  DENSITY 
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5.3.3  Effective  Body  Analysis 

The  concept  of  an  effective  body  to  replace  the  separated  region  or  the 
displacement  thickness  of  boundary  layer  is  examined  for  Model  4.  The  effective 
body  is  defined  as  a  new  body  for  which  the  calculated  inviscid  surface  pressure 
agrees  with  the  measured  data  for  the  original  body  geometry.  Because  the 
separation  bubble  is  so  large  for  Model  4  that  the  new  body  will  have  the  same 
order  of  coordinate  perturbation  in  both  x  and  y  directions,  therefore  a  further 
assumption  that  the  pressure  is  constant  along  the  normal  direction  of  the 
original  body  (a  boundary  layer  like  assumption)  within  the  viscous  layer  is  made. 
The  new  body  is  then  obtained  by  trial  and  error  until  the  surface  pressure  matches 
with  the  experimental  data  as  shown  in  Fig.  11a.  In  Fig.  lib,  a  comparison  is 

made  between  the  effective  body  shape  and  the  edge  of  separated  region  obtained 
from  the  flowfield  picture  (Fig.  11a).  Also  included  in  the  figure  is  the 
measured  shock  location  and  the  computed  ones  using  both  the  actual  and 
effective  body  shapes.  Good  agreement  is  obtained  between  the  effective  body 
calculation  and  measurements.  The  implication  of  the  effective  body  analysis  is 
that  the  separated  region  may  be  considered  as  a  solid  portion  of  the  body  and 
inviscid  solution  for  the  effective  body  give  a  better  flowfield  than  does  the 
actual  body. 
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5.3.4  Mach  Number  Effects 

In  actual  flight,  the  freestream  Mach  number  is  higher  than  what  have  been 
calculated.  To  see  the  Mach  number  effects  on  the  surface  pressure  distribution, 
the  H»  ■  14  case  is  compared  to  that  of  Moo  =  5  for  Model  4  as  shown  in  Fig.  12 
and  no  significant  change  is  shown.  The  procedure  used  to  obtain  M»  =  14  results 
is  by  continuously  increasing  the  freestream  Mach  number  (the  sequence  is 

■  6,8,10,12,14)  starting  from  the  Ha  =  5  solution  obtained  previously  without 
changing  the  number  of  grid  points. 


FIGURE  12.  CALCULATED  SURFACE  PRESSURE  COEFFICIENT  AT  DIFFERENT  MACH  NUMBER 
5.4  Coupling  with  Afterbody  Code 

A  coupling  of  the  inviscid  nosetip  calculation  with  an  existing  NSWC 
afterbody  codell  has  been  accomplished.  An  initial  plane  flowfield  data  is  stored 
in  tape  with  the  appropriate  format  to  be  accepted  by  the  afterbody  code.  Two 
examples  of  this  type  calculation  are  given  in  Fig.  13  and  14.  In  Fig.  13,  the 
surface  pressure  distribution  for  a  blunted-nose-cone  (LAM)  at  Moo  =  9.8  is 
compared  to  the  experimental  data.  The  initial  plane  is  located  at  x/RN  =  1.3 
and  calculation  covers  a  body  length  of  30.  The  agreement  with  the  experimental 
data  is  satisfactory.  In  the  same  figure,  a  sphere-cone  result  is  also  plotted, 
this  shows  the  influence  of  nose  shape  on  the  afterbody  pressure  distribution. 

In  Fig.  14,  the  afterbody  pressure  distributions  for  model  4  at  H»  =  5  are 
plotted  for  two  different  initial  plane  flowfields  obtained  previously;  i.e.  the 
real  body  and  the  effective  body.  The  pressure  and  the  axial  velocity  at  the 
initial  plane  x  =  4.32  in.  are  also  shown  in  Fig.  14.  As  shown  in  Fig.  14,  the 
influence  of  initial  plane  flowfield  on  the  afterbody  pressure  distribution  is 
limited  to  about  4  in.  downstream.  The  agreement  with  experimental  data  is  not  as 
good  as  the  previous  example.  The  obvious  reason  is  that  the  nosetip  flowfield 
is  not  well  calculated  because  of  flow  separation. 
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FIGURE  13.  COMPARISON  OF  INVISCID  SOLUTION  AND  EXPERIMENT  FOR  AFTERBODY  SURFACE  PRESSURE 
OF  LAM  NOSETIP  CONE  AT  M  -  9.8 


SOLUTION  AT  X=4.32 


FIGURE  14.  COMPARISON  OF  INVISCID  SOLUTION  AND  EXPERIMENT  FOR  AFTERBODY  SURFACE  PRESSURE 
OF  A  CONE  WITH  MODEL  4  NOSETIP  AT  M  =5 
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CHAPTER  6 

PROBLEMS  IN  VISCOUS  FLOW  CALCULATION  USING  KCL  CODE 


6.1  Sphere-cone 

To  verify  the  code  for  viscous  flow  calculation  based  on  the  thin-layer 
approximation  of  the  Navier-Stokes  equations  (or  thin  layer  theory) ,  laminar  flow 
over  a  sphere-cone  at  ^  =  5.92,  Reoo  =  10®  and  Tw/Too  =  4.4  was  computed.  The 
grid  used  was  32(J)  x  28(K)  with  8  =  1.005  (Eq.  4.2).  The  steady  state  solution 
was  obtained  in  400  time  steps  with  a  Courant  number  of  75  (the  non-dimensional 
shock  speed  is  in  the  order  of  10“^) .  As  shown  in  Fig.  15  the  calculated  results 
for  heat  transfer  in  term  of  Stanton  number  over  the  surface  is  compared  to  the 
measured  data  reported  in  Ref.  7.  Also  plotted  in  Fig.  15  is  the  boundary  layer 
calculation  using  Cebeci-Smith 's  boundary  layer  codeias  given  in  Ref.  7. 

It  is  seen  that  the  K-C-L  code  gives  higher  ST  value  up  to  30%  than  predicted 
by  the  boundary  layer  theory.  It  should  be  pointed  out  that  for  the  hemisphere- 
cone  case,  the  flow  is  fully  attached  and  the  surface  pressure  agrees  well  between 
the  inviscid  solution  (obtained  from  an  existing  NSWC  code^Z  and  was  used  in  the 
boundary  layer  calculation)  and  the  laminar  solution  given  by  the  K-C-L  code. 
Therefore,  one  would  expect  good  agreement  in  heat  transfer  results.  As  described 
in  Section  7,  after  a  major  modification  of  the  K-C-L  code  by  rewriting  all  the 
viscous  subroutines  according  to  the  analysis  presented  in  Sections  2  and  3,  the 
new  results  indeed  agree  well  with  the  boundary  layer  calculation. 

An  algebraic  turbulence  model  developed  by  Baldwin  and  Lomax^  was  incorporated 
in  the  original  code.  To  test  out  the  turbulent  solution,  the  final  laminar 
solution  with  Eg  =  0.1  was  used  as  the  initial  condition.  It  was  found  that  the 
turbulent  calculation  would  not  converge  if  the  implicit  discipation  terms  are 
set  to  zero,  i.e.,  El  =  0  (Note  that  although  Ref.  1  mentioned  about  the  imposing 
of  implicit  dissipation  terms,  but  these  terms  were  not  appeared  in  the  code 
received).  Only  after  the  implicit  dissipation  terms  were  added  into  the  code 
with  eI  =  3eE,  the  shock  speed  converges  well.  For  hemisphere-cone,  a  turbulent 
solution  was  obtained  in  200  time  steps  with  the  non-dimensional  shock  speed  in 
the  order  of  10”3.  As  shown  in  Fig.  15,  the  Stanton  number  increases  significantly 
for  turbulent  flow  as  compared  to  the  laminar  solution.  No  measured  data  is 
available  for  comparison.  It  should  be  pointed  out  that  the  values  of  surface 
pressure  obtained  from  the  laminar  and  turbulent  solution  agree  to  two  digits. 

A  minor  modification  in  the  distribution  of  eE  values  has  been  added  into  the 
code.  Instead  of  a  uniform  distribution  of  eE  over  all  the  grid  points,  it  is 
linearly  reduced  from  eE  at  the  shock  (k  ■  K)  to  O.leg  at  the  body  (k  =  1).  This 
is  done  for  viscous  flow  calculation  only  and  will  help  to  show  the  real  viscous 
effects  because:  (i)  the  true  viscous  terms  are  important  in  the  area  near  to  the 
wall  and  (ii)  the  flow  is  essentially  inviscid  at  the  shock  where  more  dissipation 
is  needed  to  smooth  out  the  oscillations  of  the  flow  variables  there.  When  eE  is 
linearly  reduced,  it  is  denoted  by  Eg  to  distinguish  from  the  uniform  one.  For 
hemisphere-cone,  the  effects  of  Eg  are  seen  to  be  insignificant  on  the  heat  transfer 
results  as  shown  in  Fig.  15. 
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FIGURE  15.  COMPARISON  OF  HEAT  TRANSFER  FOR  A  SPHERE-CONE  AMONG  K-C-L  CODE  SOLUTION,  MEASURED 
DATA  AND  B.L.  CALCULATION 
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6.2  Indented  Nosetips 

The  calculation  procedure  used  for  inviscid  flow  over  indented  nosetips  as 
described  in  section  5.0  was  also  applied  to  calculate  viscous  flow  over  indented 
nosetips.  Model  1  and  4  were  chosen  for  this  investigation. 

Laminar  flow  over  Model  4  was  first  calculated.  The  calculation  started  with 
a  grid  of  24  x  32  (CN  =  150,  eE  =  0.4,  eI  =  0)  for  400  time  steps  to  obtain  a 
laminar  solution  over  a  sphere  at  H»  =  5.0,  Re  =  8  x  10^/FT  and  Tw/Too  =  5.4.  The 
sphere  was  then  deformed  to  the  shape  of  Model  4  in  1800  time  steps.  The  grids 
were  then  increased  to:  (A)  58  x  32  and  (B)  56  x  48  in  another  1600  time  steps 
each  with  the  final  values  of  CN  =  50,  eE  =  0.1  for  (A)  and  eE'  =  0.3  for  (B) . 

The  calculated  surface  pressure  and  shock  locations  from  these  two  solutions  are 
close  as  shown  in  Fig.  16.  This  provides  a  self  verification  of  the  results. 

Since  grid  (B)  contains  more  points  in  the  n  direction,  its  solution  is  used  for 
comparison  with  the  measured  data  as  shown  in  Fig.  17  and  18. 

Unlike  the  hemisphere-cone,  the  turbulent  calculations  for  Model  4  encountered 
serious  difficulties.  Large  amplitude  oscillation  of  pressure  in  the  flowfield 
quickly  interrupted  the  computation.  The  value  of  CN  was  gradually  reduced  and  the 
value  of  eI  was  increased  (eE'  =0.3  was  maintained).  At  CN  =  2  and  eI  =  6,  it 
was  possible  to  run  for  200  time  steps  with  the  non-dimensional  shock  speed 
converging  to  a  value  of  0.04.  The  shock  speed  then  starts  to  increase  slowly 
but  steadily.  A  further  increase  of  eI  up  to  12  would  not  help  to  obtain  a 
converged  solution.  Thus,  the  solution  before  the  shock  speed  started  to  increase 
is  shown  in  Fig.  17  and  18  for  comparison. 

Fig.  17  shows  the  comparison  of  shock  location  between  calculations  and 
experiments.  The  inviscid  shock  layer  is  thinner  around  the  indented  region  as 
expected.  The  laminar  and  turbulent  solutions  for  shock  location  are  almost 
coincide  and  fall  in  between  the  inviscid  solution  and  the  measured  data.  The 
primary  separation  bubble  indicated  by  the  laminar  solution  is  smaller  than 
observed  experimentally.  The  laminar  separation  point  of  the  primary  separation 
bubble  is  at  the  downstream  end  of  the  expansion  corner,  but  the  flow  picture 
(Fig.  9a)  shows  that  the  flow  separates  immediately  at  the  beginning  of  the 
expansion  corner.  In  general,  the  effects  of  turbulence  is  to  move  the  separation 
point  toward  downstream  and  the  separation  bubble  will  be  smaller.  This  fact 
suggests  that  the  discrepancy  shown  in  Fig.  17  is  not  because  of  turbulence  effects 
but  from  some  other  sources  which  are  not  correctly  simulated  in  the  numerical 
solution.  Also  the  laminar  solution  indicates  that  there  is  a  secondary  separation 
bubble  within  the  primary  separation  bubble  as  sketched  in  Fig.  17.  The  secondary 
separation  bubble  was  not  reported  in  Ref.  8. 

In  Fig.  18,  the  surface  pressure  distribution  obtained  from  the  inviscid, 
laminar  and  turbulent  solution  are  compared  to  the  measured  data.  It  is  noted 

that  the  viscous  solutions  compare  better  with  the  measured  data  than  the  inviscid 

curve.  The  region  around  the  expansion  corner  S/L  'v  0.26  -  0.7  where  the  inviscid 

and  viscous  solutions  are  seen  to  agree  well  (i.e.,  no  flow  separation)  but  are 

lower  than  the  measured  data.  The  dip  in  the  pressure  curve  in  the  region 
S/L  ^2.5  (where  the  secondary  separation  bubble  starts)  is  not  shown  by  the 
measured  data. 
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FIGURE  17.  COMPARISON  OF  SHOCK  LOCATION  BETWEEN  K-C-L  CODE  SOLUTION  AND 
EXPERIMENT  FOR  MODEL  4 


FIGURE  18.  COMPARISON  OF  SURFACE  PRESSURE  BETWEEN  K-C-L  CODE  SOLUTION  AND  EXPERIMENT  FOR  MODEL  4 
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The  most  troublesome  result  obtained  from  the  K-C-L  code  for  Model  4  is  the 
heat  transfer  rate.  For  S/L  >  0.68  (after  flow  separation)  the  heat  transfer  rate 
becomes  negative  which  is  not  realistic. 

For  Model  1,  there  is  a  sharp  expansion  corner  with  a  radius  of  0.062  inch. 

Three  grid  points  were  used  to  cover  the  corner  as  was  done  for  the  inviscid 
calculation  and  a  total  grid  points  of  33  x  32  were  used  for  the  viscous 
calculation.  Only  laminar  solution  can  be  obtained.  As  shown  in  Fig.  19  and  20, 
while  the  inviscid  solution  agrees  reasonably  well  with  the  measured  data  for  both 

the  shock  location  and  the  surface  pressure,  the  laminar  solution  is  very  poor. 

The  reason  is  that  the  flow  separates  immediately  after  the  corner  and  form  a 
large  primary  separation  bubble  as  shown  in  Fig.  19.  Within  the  primary  separation 
bubble,  there  is  also  a  secondary  separation  bubble  around  the  location  of 
compression  turn.  As  a  result  of  the  primary  separation  bubble,  the  laminar  shock 
layer  becomes  thicker  near  the  separation  bubble  and  thinner  afterwards  as  compared 
to  the  measured  data.  The  surface  pressure  obtained  from  the  laminar  solution  looks 
entirely  wrong  as  shown  in  Fig.  20. 

It  was  not  possible  to  obtain  a  turbulent  solution  for  Model  1,  not  even  one 
like  that  of  Model  4.  The  obvious  reason  is  that  the  laminar  solution  is  too  far 

off  from  the  measured  data,  which  is  assumed  to  be  close  to  the  turbulent  solution, 

therefore  the  starting  flowfield  is  too  poor  to  carry  through  the  calculation. 

With  all  the  troubles  in  simulating  the  viscous  flowfield,  particularly  the 
temperature  field,  an  effort  was  made  to  repeat  the  results  given  in  Fig.  5-8  of 
Ref.  1  for  a  hemisphere-cylinder  with  adiabatic  wall  at  M»  =  2.94  Reoo  =  2.2  x  10^ 
and  To  =  293°K.  It  was  found  not  possible  to  duplicate  the  temperature  results 
given  in  Fig.  7  of  Ref.  1  with  the  K-C-L  code,  but  the  obtained  surface  pressure 
and  shock  shape  do  repeat  well.  Therefore,  ,a  reanalysid  was  carried  out  for  the 
entire  calculation  procedure  as  shown  in  Sections  2  and  3  and  also  all  the  viscous 
subroutines  were  rewritten.  As  given  in  the  next  section,  a  significant 
improvement  of  the  results  of  temperature  field  is  obtained. 
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FIGURE  19.  COMPARISON  OF  SHOCK  LOCATION  BETWEEN  K-C-L  CODE  SOLUTION  AND 
EXPERIMENT  FOR  MODEL  I 
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FIGURE  20.  COMPARISON  OF  SURFACE  PRESSURE  BETWEEN  K-C-L  CODE  SOLUTION  AND 
EXPERIMENT  FOR  MODEL  1 
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CHAPTER  7 

VISCOUS  FLOW  CALCULATIONS  USING  THE  NEW  CODE 

A  modification  of  the  K-C-L  code  by  rewritten  all  the  viscous  subroutines  has 
been  accomplished  and  calculations  were  performed  for  laminar  flows  over  a 
hemisphere-cylinder  and  a  sphere-cone  using  the  meshes  A,  B,  C  and  D  given  in 
Table  2  (section  4.2).  Comparison  of  results  is  described  in  this  section. 

Because  of  the  termination  of  support,  calculations  for  indented  nosetips  and 
incorporation  of  turbulent  model  have  not  been  carried  out. 

7.1  Hemisphere  -  Cylinder 

The  results  for  hemisphere-cylinder  with  adiabatic  wall  at  Nko  of  2.94, 

Reoo  of  2.2  x  1()6  and  To  of  293°K  are  given  in  Figs.  21  to  23.  Figure  21  shows 
the  temperature  distribution  over  the  body  surface.  It  is  seen  that  the 

results  obtained  from  Mesh  B  and  C  agree  quite  well  but  not  that  given  by  Mesh  A 
which  provides  not  enough  points  in  the  viscous  layer.  As  shown  in  Table  2,  the 
distribution  of  mesh  point  differs  significantly  between  Mesh  B  and  C  and  the 
solutions  agree  well  (temperature  is  a  more  sensitive  variable  than  other 
variables).  Thus,  it  is  necessary  to  provide  sufficient  grid  points  to  resolve 
the  viscous  effects  near  the  surface,  such  as  that  given  by  Mesh  B  or  C.  Also 
plot  in  Fig.  21  is  the  solution  of  Viviand  and  GhazzilO  who  solved  the  full 
Navier-Stokes  equation  and  slight  differences  are  found  in  the  area  near  the 
shoulder  between  his  and  the  present  solution.  It  was  unable  to  produce  the 
results  of  Kulter  et  al  as  given  in  Fig.  7  of  Ref.  10  from  the  copy  of  computer 
code  supplied  by  him.  The  result  of  T/T<»  given  by  Kutler's  code  using  Mesh  B 
is  shown  in  dotted  line  which  is  obviously  wrong.  Kutler  et  al  also  indicated 
that  temperature  is  a  mesh  dependent  variable,  it  seems  not  so  as  seen  from  the 
present  results  of  Mesh  B  and  C  given  in  Fig.  21. 

The  temperature  profile  at  three  stations  obtained  from  Mesh  B  and  C  are  shown 
in  Fig.  22,  the  agreement  with  the  solution  of  Viviand  and  Ghazzi  at  slightly 
different  station  (given  in  parenthesis)  is  good.  The  shock  shape  and  surface 
pressure  are  plotted  in  Fig.  23.  It  is  interesting  to  note  that  these  two 
quantities  are  not  as  sensitive  as  temperature  and  all  solutions,  including  the 
one  obtained  using  the  K-C-L  code,  agree  well. 

7.2  Sphere-Cone 

The  results  for  sphere-cone  with  cone  half  angle  of  9.75  deg  and  with 
isothermal  wall  of  Tw/Too  *  4.4  at  Moo  =  5. 92  and  Reo,  =  10^  are  shown  in  Fig.  24-26. 
In  hemisphere-cylinder  calculation,  it  is  noted  that  Mesh  B  gives  as  accurate 
results  as  Mesh  C  but  with  much  less  computing  time,  thus  the  same  value  of 
B  *  1.005  is  chosen  for  sphere-cone  calculation.  In  order  to  compare  heat  transfer 
with  boundary  layer  calculation,  it  is  important  that  the  surface  pressure  used 
in  the  boundary  layer  calculation  must  be  consistent  with  that  obtained  from  the 
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FIGURE  25.  HEAT  TRANSFER  DISTRIBUTION  IN  TERMS  OF  STANTON  NUMBER  OF  SPHERE 
CONE  WITH  ISOTHERMAL  WALL  OF  Tw/T«,  =  4. 4  AT  H»  ■  5.92  AND 
Re*  -  10° 
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FIGURE  26.  COMPARISON  OF  TEMPERATURE  PROFILE  NORMAL  TO  SURFACE  BETWEEN  PRESENT 
SOLUTION  AND  BOUNDARY  LAYER  CALCULATION  FOR  SPHERE-CONE  WITH 
ISOTHERMAL  WALL  OF  Tw/T«o  *  4.4  AT  Moo  =  5.92  AND  Reco  =  106 


present  solution.  In  Fig.  24  a  comparison  is  made  for  the  surface  pressure 
between  the  present  solution  and  the  inviscid  surface  pressure  solution  (using 
a  blunt  body  code  developed  at  NSWC)2*  used  in  the  boundary  layer  calculation. 

The  agreement  is  excellent  except  at  the  shoulder  where  slight  difference  is 
shown.  The  Stanton  number  distribution  is  given  in  Fig.  25.  It  is  seen  that  the 
agreement  between  present  solution  and  the  boundary  layer  calculation  is 
surprisingly  good  up  to  the  shoulder,  from  there  on  slight  difference  is  shown. 
The  solution  obtained  with  K-C-L  code  gives  overall  higher  values  of  ST.  Because 
of  wide  spread  in  the  experimental  data,  all  calculations  seem  to  fall  within  the 
experimental  error  band. 

A  comparison  of  temperature  profile  at  several  stations  between  present 
calculation  and  that  of  boundary  layer  calculation  (Station  given  in  parenthesis) 
is  given  in  Fig.  26.  The  agreement  is  good.  At  station  S/Rn  =  1.4,  i.e. 
the  shoulder,  again  the  difference  is  larger.  The  good  agreement  of  present 
calculation  with  boundary  layer  calculation  for  the  temperature  field  is  a  good 
check  of  the  code  since  the  surface  pressure  distributions  in  both  cases  agree 
well. 


21Courant,  R. ,  and  Hilbert,  D.,  "Methods  of  Mathematical  Physics,"  Vol.  2, 
Interscience,  New  York,  1062,  pp.  558-605. 
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CHAPTER  8 

SUMMARIES  AND  RECOMMENDATIONS 

Extensive  applications  of  the  K-C-L  code,  which  solves  the  Euler  or 
Navier-Stokes  equations  with  thin-layer  approximation,  have  been  performed 
for  inviscid  and  viscous  flows  over  smooth  (hemisphere-cylinder  and 
sphere-cone)  and  indented  (given  in  Section  4.1  or  Fig.  3)  nosetips  at 
hypersonic  speed.  Comparisons  of  calculated  results  and  measured  data  are 
made  for  surface  pressure,  shock  location  and  heat  transfer  rate. 

The  inviscid  portion  of  the  K-C-L  code  works  well  for  smooth  nosetips.  For 
indented  nosetlp  calculations,  a  special  calculation  procedure  has  been 
developed  in  order  to  run  the  code  and  reasonable  solutions  for  surface 
pressure  are  obtained  as  compared  to  the  measured  data  when  the  separation 
bubble  is  small. 

The  calculated  inviscid  flowfields  for  the  series  of  indented  nosetips  under 
investigation  Indicate  that  there  is  no  strong  embedded  shock  or  slip 
surface  within  the  shock  layer  for  the  cases  investigated.  The  most  serious 
aerodynamic  problem  for  this  series  of  nosetips  is  flow  separation,  which 
requires  a  solution  to  the  Navier-Stokes  equations. 

Because  of  the  poor  results  of  the  temperature  field  given  by  the  K-C-L 
code,  a  reanalysis  of  the  complete  calculation  procedure  has  been  carried 
out  as  given  in  Sections  2  and  3  and  all  the  viscous  subroutines  have  been 
rewritten. 

The  modified  code  gives  good  results  for  the  temperature  field  as  demonstrated 
for  the  cases  of  laminar  flow  over  hemisphere-cylinder  and  sphere-cone.  It 
is  believed  that  the  modified  code  can  provide  reliable  prediction  of  laminar 
flowfield  for  nosetips  without  flow  separation. 

In  order  to  predict  indented  nosetip  with  large  separation  bubble,  it  is 
necessary  to  extend  the  present  code  to  include:  (i)  a  solution  to  the 
full  Navier-Stokes  equations  with  considerable  more  grids  available  to  cover 
the  separated  region  and  (ii)  a  good  trubulence  model  (to  be  searched  and 
tested)  since  the  flow  is  likely  to  be  turbulent  after  separation. 
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APPENDIX  A 

JACOBIAN  MATRICES  FOR  A,  B,  K,  L,  M,  N,  P,  Q 
The  Jacobian  matrices  used  in  Eq.  2.7  are  listed  in  the  following: 


kicr 

uCkjU+^v) 

kjj-kj  (y~2)u  + 
kLu+k2v 

k2cr 

v(klU+k2v) 

kfV-kj  (Y~0  u 

-kj  (Y-l)v 


+k2u 


kD-(Y-2)k2v 

+k1u+k2v 


0 


kj (Y-1) 


k2(Y-l) 


?  'ci)kr 

Y“l)  (kjU+k2v)  u 
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L  or  M  ■  _i 
JP 


0 

0 

0 

0  ' 

v(  ) 

X 

0 

-(  >x 

0 

v(  )y 

0 

-<  >y 

0 

2v[u(  )x+v(  )y] 

-v(  )x 

-[u(  )X+2v(  )y] 

0 

(A.  3) 


where  (  )  =  C  for  L  and  (  )  =  r)  for  M. 


P  or  0  -- 


0 

0 

0 

-(  )2u-(  )3V 

(  )2 

<  >3 

> 

W 

1 

o 

cn 

/-s 

(  >3 

<  >4 

-(  >,*-**> 

-(  )5u+(  )2u 

-(  >5V+(  >4v 

-(  )2u2-2 (  )3uv 

+(  )3V 

+(  )3u 

-(  )*v2 

(A.4) 


where  (  )  =  b  for  P  and  (  )  =  c  for  Q. 


0 

0 

0 

0 

1 

0 

0 

0 

0 

jpy 

3c6v 

0 

-3c6 

0 

2c6v2 

Hi 

T*2VT1 

-|d2un-2c6v 

0 

(A.  5) 


A- 2 


wlu, 
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APPENDIX  B 

DERIVATION  OF  EQ.  (3.9) 


To  obtain  Eq.  (3.9)  through  characteristic  analysis  is  briefly  described  in 
this  appendix.  The  background  of  characteristic  theory  and  characteristic 
compatibility  conditions  may  be  found  in  Ref.B-1  and  their  applications  in  fluid 
dynamic  may  be  found  in  Ref.  B-2. 


Equation  (1.1)  for  inviscid  flow  is 


Q.  +  F  +  F  +  F--A  =  0 

ut  .x  ry  y 


(B.l 


It  is  preferred  to  change  the  dependent  variable  from  y  to  Q, 


Q  = 


/  P 

u 
v 
P 


(B.  2 


The  resulting  equation  is 


®t  +  B05«  +  CoQj,  +  D0  -  0 


(B.3 


where 


r 

p 

0 

0 

0 

u 

0 

1/P 

0 

0 

u 

0 

\° 

YP 

0 

u 

/v 

0 

P 

0 

1 0 

V 

0 

0 

0 

0 

V 

1/P 

0  Yp  v 
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The  characteristic  matrix  for  Eq.  (3.3)  is 


A0*  (^i.  X2,  V  =  X1T  +  X2Bq  +  X3Cc 


X2p  X3P  0 
a0  0  X2/p 

0  a0  3/p 

X2YP  X3yp 


where  aQ  =  Xj  +  uX2  +  VX3 

With  the  transformation  of  Eq.  (1.4),  Eq.  (B.3)  becomes 

+  B1Q£  +  ciQn  =  ~Do 

where  B^  =  +  +  £y  Co 


Ci  =  V  +  r^0  +  ny  Cc 


The  characteristic  matrix  for  Eq.  (B.5)  is 


A^  (Xp  X2,  X3)  =  X-^I  X2®1  X3(-l 


(B.  4) 


(B.4a) 


(B.5) 


=  V  +  x2A*  (Ct,  Cx,  Sy)  +  VoK’VV 

=  A*^  +  ?tX2  +  ntX3,SxX2  +  \X3,CyX2+riyX3) 


Hence 


-  A0  (X^  *  X2 ,  X3) . 


xi  -  \  +  ct’x2  +  nt) 
X2  =  ^xX2  +  nxX3 


(B.  6) 


(B.  7) 


X3  *  ?yX2  +  nyX3 
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The  characteristic  condition  is 


det.  A*  =  aQ2  [aQ2  -  a2(X22  +  X32)]  =  0 

(B.8 

0 

» 

t>° 

(B.8 

cr  =  +  aA02  +  X,2 
o—2  3 

(B.8 

The  characteristic  curves  <j>(x,  C.  n)  corresponding  to  the  four  distinct 
characteristic  conditions  of  Eq.  (B.8)  have  slopes: 


dE 


S  !  2  =  ct  +  ^  +  ^y 


(B.9a 


dn 

dx  ,  „  =  nt  +  urix  +  vny 


1,2 


(B.9b 


for  a  =0  and 
O 


(fjf)  =  h  +  u?x  +  v?y  -  [<SX2  +  ^y2)  *2  +  <*x\  + WV 

\  /3  j  A  o 


(B.1C 


(s)3  4  ■  \  +  unx +  '"V  -  ^  +  W^2  +  +  \2)V  <B-1C 


A  A 

for  Oq  =  +  aA^  +  X^  .  On  the  boundary  n  =  0,  the  physical  boundary  condition 
requires  v  =  nt  +  ur)x  +  vr)y  =  0,  the  slopes  of  (dr|/dx)3  2  are  zero  and  that 


W\  _ 

(dT/3,4  =  +  a2yTFT 


==  [nx5x  +  TL?  )X2  +  (nx2  +  Oy2)X3l.  Hence,  the 

2“  ’  X3 


admissible  characteristic  conditions  are  1,  2  and  3  (see  sketch) 

*1 


Sketch  of  Bicharacteristic 
Directions 


,<X0  =  0 
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To  obtain  the  compatability  conditions  corresponding  to  the  three  admissible 
characteristic  conditions  1,  2  and  3,  one  finds  the  left  null  vector  lt^,  x2  and 
by  requiring  J  "  J 


Then, 


£  *  A1  =  0 

=  (0»  »  —  ^3^x»0) 

S.2  =  (  P»  0,  0,  p) 


(B.ll) 

(B. 12a) 
(B. 12b) 

HyY?  ,1  (B. 12c) 

a/lx2  +  ny2 


where  the  case  of  A2  =  0  and  4  0  is  considered.  The  compatability  condition  is 
then  obtained  by 

l  •  Aj  (1,0,0)  Qt  +  %  •  A*  (0,1,0)  Qx  +  I  *  A*(0,0,1)  Qy  +  1  •  Dq  =  0  (B.13) 

It  is  noted  that  the  compatibility  conditions  for  and  &2  are  the  equations  for 
a  fluid  particle  or  streamlines, 

For  Jtj: 

ny  (ut  +  uux  +  vuy  +  —  )  +  nx(^r  +  uvx  +  wy  +  lx.  =  0  (B.  14) 

which  is  automatically  satisfied  when  the  momentum  equations  are  satisfied,  and 


Pt  +  (Pu>x  +  (pv)y  _  ~2<Pt  +  UPX  +  vpy)  =  0  (B.  15) 

which  is  the  principle  of  constancy  of  particle  entropy.  Thus  Eqs.  (B.14)  and 
(B.15)  describe  a  particle  path.  The  last  compatibility  condition  is 


pa 


p_  = - >  r.  i  ■  ...  (n  +  up  +  vri  ) 

T  A  2  ~  ~  2  «  XT  V 

x  y 


+  aA]  2  +  r)  2  p_  -  pa^(£  ur  +  £  v_  +  r)  u  +  riv) 
x  y  rn  x  t,  y  t,  x  n  y  n 


v 

y 


-  up^  +  -j  == — - [nx(puu^  +  ^xp^)  +  ny(puv^  +  Cyp^)l  (b*  16) 

Therefore,  the  pressure  on  body  surface  can  be  integrated.  Once  the  pressure  is 
obtained  the  rest  of  the  flow  variables  may  be  determined  by  the  isentropic 
relations  as  described  in  section  3, 
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APPENDIX  C 

CALCULATION  OF  6  FOR  INDENTED  NOSETIP  SHAPES 

For  indented  nosetip  shapes  given  in  Fig.  3,  expressions  for  the  body  point 
(x^(0),  y^(0))  and  distance  of  deformation  6(0)  are  listed  in  the  following. 

Given:  control  points  (XPn,  YPn) ,  n  *  1  to  7 

Rn:  n  =  1  to  4 

Bn:  n  =  1  to  3 

The  centers  of  circular  arc  for  expansion  comer,  compression  turn  and 
expansion  shoulder  are: 

Xq-^  =  R^  (for  flat  nose  R-^  -*■  <») 

Xq2  =  XP 2  +  R2  *cos  Si 

Xq2  =  Y?2  -  R2  *sin  0^ 

Xq3  =  XP3  -  R3  *cos  0^ 

Y03  =  YP3  +  R3  sin 

Xq4  =  TCP ^  +  R4  cos  02 
Y04  =  ^*5  “  ^4  8^-n  02 

Let  the  center  of  the  sphere-cone  be  located  at  (XQO,  0).  The 

value  of  Xqq  (radius  for  the  sphere)  must  be  greater  than  XPg.  Then  the  0n  for 
each  control  point  are  given  by: 

0X  =  tan'1(YP1/(X0i  -  XP^) 

02  =  tan'1(YP2/(X01  -  XP2)) 

03  -  tan"1(YP3/(X01  -  XP3)) 

04  =  tan_1(YP4/(X01  -  XP4)) 

05  =  tan"1(YP5/(X01  -  XP5)) 

06  =*  tan‘1(YP6/(X01  -  XPg)) 

C-l 
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e2  <  9  1  63! 


03  <  6  <  64: 


=  (XP2  +  A  tan  81)/(1  +  tan  8  tan  8^) 


A  =  tan  8  -  YP0 
00  2 


B  +  A1  -  C 
A  =  XQO  tan  6  -  Y03 
B  =  (X03  +  A  tan  0)  cos20 
C  =  (Xfto2  +  A2  -  Ro2)  cos20 


04  <  9  1  05:  Xfi  =  03*4  +  A  tan  62) /(I  +  tan  6  tan  62> 

A  =  X  tan  0  -  YP, 

00  4 


05  <  9  <  06: 


Xg  =  B  -  /B2-C 

A  =  Xoo  tan  9"Y04 
B  =  (Xq4  +  A  tan  0)  cos20 

C  =  (Xq42  +  A2  -  R42)  cos20 


)&  <  0  <  0m:  XB  =  (XP6  +  A/tan  B3)/(l  +  tan  0/tan  83) 


A  =  X_„  tan  0  -  YP, 


The  corresponding  Y^  and  6  are  given  by: 


Yb  =  (Xoo  '  x)  tan  0 


6  =  Xqq  [1  -  (I-7-  ) / cos  0] 


Therefore,  a  sphere-cone  with  nose  radius  equals  to  Xqq  and  cone  half  angle 
9C  =  63  use^  as  t^ie  initialized  shape  and  let  it  deformed  to  the  desired 
indented  nosetip  shape. 


C-3 
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APPENDIX  D 

PROGRAM  DESCRIPTION  AND  OPERATING  MANUAL 

This  computer  program  contains  one  main  program  called  NOSETIP  and  twenty  four 
subroutines.  The  main  program  NOSETIP  serves  as  a  flow  chart  as  given  in  Fig.  D1 
which  also  describe  the  structure  of  the  computer  program.  A  brief  description 
of  the  function  of  each  subroutine  is  given  in  Table  D2.  A  listing  of  important 
Fortran  symbols  is  given  in  Table  D3  and  the  complete  listing  of  the  program 
follows. 

The  operating  manual  is  best  described  by  the  input  data  cards  for  different 
cases  that  the  computer  program  can  handle.  This  is  given  in  Table  D4.  The 
corresponding  examples,  seven  all  together,  and  the  output  data  are  separately 
described  immediately  after  Table  D4. 
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Three  different  intializations : 

(1)  Fresh  start;  (2)  Restart;  and 
(3)  Adding  mesh  points.  Called 
subroutines:  CTXY,  ETATB,  GRID, 
JACOB,  SHAPE,  INTERP.  Print  out 
input  data  and  free  stream 
conditions. 

Prints  out  starting  flowfield. 

Determines  starting  step  size 

Computes  starting  right  hand  side 
of  Eq.  (2.6). 

Print  out  starting  residue 
information. 

Advances  shook  points  explicitly. 

Advances  body  points  explicitly. 
Called  subroutine  TRIB. 

Integrates  equations.  Called 
subroutines:  RHS ,  RESIDU,  DISSIP, 
EFCON ,  ABMATX ,  LBLTRA,  LBLTRB, 
LUDEC,  BTRI ,  VSMATB,  VSRHSB 

CHANGES  MESH  AND  COMPUTES  new 
metric  terras. 

Computes  new  step  size  every  25 
steps. 


Prints  out  solutions  or  stores  data 

(1)  Flowfield  information 

(2)  E  and  F  conservative  variables 

(3)  store  solution  for  restart  on 
tape  2 

(4)  store  solution  for  afterbody 
calculation  on  tape  3. 

(5)  Residue  information 

(6)  Heat  transfer  information 


FIGURE  Dl.  FLOW  CHART 
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Table  D2.  Brief  Description  of  Subroutines 

Subroutine  Function 

Name 


ABMATX(J,K, I) : 

Computes  coefficients  for  matrix  A(I=1)  and  B(I*2)  as 
given  by  Eq.  (A.l)  at  nodal  point  (J,K) 

BNDRY : 

Advances  explicitly  the  body  boundary  condition  in  time 
see  body  points  in  section  3.0. 

BTRI(IL.IU): 

Solves  a  block  tridiagonal  system  of  equations.  IL  and 

IU  denote  the  starting  and  finishing  indices. 

CTXY (XAA,  YAA, 

DED,  BF,  CT): 

Obtains  (x,y)  6  and  bg  values  corresponding  to  a 
given  9  (=CT  in  radian)  value  for  the  nosetip  shape 
described  in  Appendix  C. 

DISSIP: 

Computes  the  explicit  dissipation  terms. 

EFCON(J,K, I) : 

Computes  the  conservative  flux  terms  given  by  E  (with  1*1) 
and  F(with  1=2)  at  nodal  point  (J,K). 

EIGEN: 

Computes  step  size  for  given  CN  by  Eq.  (4.5). 

ETATB (ET ,  CF , KMAX) : 

Computes  the  mesh  distribution  in  r)  direction  according 
to  Eq.  (4.2a)  with  CF  =  B. 

GRID: 

Sets  up  a  grid  in  the  physical  plane  at  the  starting 
of  a  new  computation. 

INITIA: 

Initializes  the  flowfield  and  prints  out  the  freestream 
conditions  and  inputs. 

INTEGR : 

Integrates  the  system  of  equations  in  time.  First 
sweep  over  £  direction  and  then  n  direction. 

INTERP : 

Interpolates  the  flow  variables  when  more  grid  points 
are  needed  at  the  beginning  of  a  continuation  run. 

JACOB : 

Computes  the  matrices  given  by  Eq.  (2.6). 

LBLTRA(K) : 

Computes  all  the  block  matrices  for  all  J  arrays  in 

E,  sweep,  see  Eq.  (2.11). 

LBLTRB(J) : 

Computes  all  the  block  matrices  for  all  K  arrays  in 

H  sweep,  see  Eq.  (2.12). 

LUDEC(A) : 

Computes  L-U  decomposition  elements,  for  2D  problem 

A  is  (4  x  4)  matrix. 

D-3 
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OUTPUT (L) : 

Prints  out  results  and  other  information  identified  by 

L  (see  operating,  manual) . 

RESIDU: 

Computes  residues  at  each  nodal  point. 

RHS: 

Computes  the  terms  in  the  right  hand  side  of  Eq.  (2.6) 
without  the  viscous  part. 

SHAPE: 

Reads  and  writes  the  control  parameters  for  the  nosetip. 

SHOCK: 

Advances  explicitly  the  shock  boundary  condition  in 
time,  see  shock  point  in  section  3.0. 

TRIP: 

Solves  tridiagonal  system  of  equations. 

VSMATB (J) : 

Computes  all  the  viscous  terms  in  the  block  matrices 
for  all  k  array  in  n  sweep,  see  Eq.  (2.12). 

VSRHSB: 

Computes  all  the  viscous  terms  in  the  right  hand  side 
of  Eq.  (2.6) 

EM 
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Table  D3.  Listing  of  Important  Fortran  Symbols 


Fortran 

Symbol 


A(J,4,4) 

Matrices  for 
of  (2.12). 

AB(4 ,4) 

Elements  for 

B(J ,4,4) 

Matrices  for 

C(J,4,4) 

Matrices  for 
(2.12). 

Cl 

Description 

or  ^  in  Eq.  (2.11) 

matrices  A  or  B,  Eq.  (A.l) 
Tj  or  Tj  in  Eqs.  (2.11)  or 

ipj+1  or  in  Eqs.  (2.11) 


Principal 

Defining 

Routine 

Common 

Block 

LBLTRA 

LBLTRB 

COM4 

ABMATX 

COM3 

(2.12) 

LBLTRA 

LBL'i  KB 

COM4 

of 

LBLTRA 

LBLTRB 

COM4 

VSMATB 

1 

VISC 

i 

C2 

i(At/Re)C3/J 

C3 

i(At/Re)C4/J 

C4 

i(At/Re)C5/J 

C5 

^At/Re)-!  . 

C6 

I  ’  ^(At/Re)  \  ny 

C7 

\  ~(At/Re)  i  nx 

CC 

CF 

CINF 

CMUKAP 


CN 


x  stretching  parameter  in  cone  portion 

GRID 

C0M1 

8  in  Eq.  (4.2a) 

ETATB 

C0M1 

free  stream  sound  speed 

INITIA 

C0M1 

(\i/\ ioo)  =  +  CVIS)/(r|-  +  CVIS) 

Aco  ^00 

VSMATB 

VI SK 

Courant  number,  Eq.  (4.5) 

INITIA 

C0M1 

CS1 


1,  At  . 
y  2Re)l1 


<^j>  \ v 


VSMATB  VISC 


D-5 
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k 


► 


CS3 

CS4 

!«'/'«>  <SrljHv  +  v) 

CS5 

-J-(At/Re)(l/py2J)nx 

CS6 

i(At/Re)(l/py2J)ny 

CS7 

Y/l  At.  1 

?r  2  ReJ^  ny  py 

\ 

/ 

CVIS 

11Q/T00(°K) 

INITIA 

CVIS1 

CVIS  +  1. 

INITIA 

D(80,80) 

J 

JACOB 

DT 

AT 

EIGEN 

DETL(80) 

Deformation  5(0)  *FACTB 

FACOB 

DETT(80) 

Deformation  6(0) 

GRID 

El  INF 

Pco/Poo 

Y-l 

INITIA 

ENT 

Y 

Entropy  behind  shock  front,  (pt/p^) 

INITIA 

ET 

a(k),  Eq.  (4.2a) 

ETATB 

ETING 

Free  stream  internal  energy  of  the  gas 

INITIA 

FACTB 

Fraction  of  deformation  for  current  run 

INITIA 

FACTT 

Fraction  of  total  deformation 

INITIA 

GAM 

Y 

INITIA 

GAMM1 

(Y-l) 

INITIA 

GAMP1 

(Y+l) 

INITIA 

GAM1I 

a/Y 

INITIA 

D-6 


COMl 

com 

COM2 

COMl 

com 

COMl 

COMl 

COMl 

COMl 

COMl 

COMl 

COMl 

COMl 

COMl 

com 

COMl 
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H 

HTINF 

IAFBD 

IGEM 


IPRT 

IR1 

IT 

ITER 

ITF 

ITRAN 

ITS 

ITWA 

IVIS 

IW2 

J 

JM 


At/2 

EIGEN 

C0M1 

Free  stream  total 

enthalpy  of  the  gas 

INITIA 

C0M1 

=1,  store  data  for  afterbody  calculation 

=0,  do  nothing 

INITIA 

C0M1 

=0  uniform  points 

on  sphere  for  sphere  cone 

GRID 

1 

C0M1 

1 

=1  read  in  XB,  YB,  XS,  YS 

=2  read  in  PH(J)  and  DETT(J)  for  arbitrary 
body  shape 

=3  uniform  spacing  for  TH(J) ,  calculate 

DETT(J)  and  determine  XB,  YB  for  indented 
nosetip. 

=4  read  in  TH(J) ,  calculate  DETT(J)  and 
XB,  YB  for  indented  nosetip. 

=1  detailed  printout  from  EIGEN  INITIA  C0M1 

=0  do  nothing 

=1  read  starting  flowfield  from  tape  1  INITIA  C0M1 

=0  do  nothing 

time  step  for  current  run  INITIA  C0M1 


total  time  steps  for  current  run 

=1  printout  heat  transfer  information 
=0  do  nothing 

No.  of  stations  at  cone  portion  for  0 
to  go  to  j 

accumulated  time  steps  from  fresh  start 

=  1  Isothermal  wall 
=  0  Adiabatic  wall 

=  1  Viscous  calculation 
=  0  Inviscid  calculation 

=  1  Write  results  on  tape  2 
=  0  do  nothing. 

Index  for  £  direction 

JMAX-1 


INITIA  C0M1 

INITIA  corn 

INITIA  C0M1 


INITIA 

C0M1 

INITIA 

C0M1 
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JMAX 

JNM 

JWRIT 

K 

KM 

KMAX 

KRES 

LIP 

OMEGA 

PINF 

PRD 

PRT 

PT 

PTORT 

Q(80,80,4) 

QINF 

REY 

REYIN 

RINF 

RR(80) 

S(80,80,4) 

SINF 

SMU 

SMUIMP 


Maximum  J  INITIA 

Value  of  J  at  junction  of  sphere  and  cone.  INITIA 

Station  with  vertical  flowfield  data  to  be  OUTPUT 

stored  for  afterbody  calculation  using 

SWINT. 

Index  for  ri  direction  INTIA 


KMAX-1 
Maximum  K 

Interval  for  printout  residues  in  K-plane 

Number  of  time  steps  to  complete  the  body 
change  in  current  run. 

Body  radius;  when  IGEOM  =  3  or  4,  this 
value  is  recalculated  by  subroutine  shape. 

=  0  for  adding  mesh  points.  ^ 

Free  stream  pressure  p^,  INITIA 

Prandtl  number 

Turbulent  Prandtl  number 

Pt,  total  pressure,  Eq.  (2.22) 


Pt/pt,  Eq.  (2.22) 

1 

Vector  U 

INITIA 

Free  stream  velocity  of  <*> 

INITIA 

Reynolds  number,  Re^  (~^^~) 

INITIA 

Read  in  R^  =  poo  qoo  L/Uoo 

INITIA 

Free  stream  density, 

INITIA 

1/P 

VSMATB 

Intermediate  values  of  U  vector 

RHS 

Free  stream  entropy 

INITIA 

Coefficient  of  explicit  dissipation, 

DISSIP 

Coefficient  of  implicit  dissipation,  eT 

LBLTRA 
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TAU 

Accumulated  time,  E  At 

INTEGR 

C0M1 

TC(80) 

£  +  H  (u2  +  v2) 

P 

VSRHSB 

TH(80) 

0  array 

INITIA 

COM1 

TM 

Cone  half  angle  in  deg 

INITIA 

COM1 

1 

i  V 


X(80,80) 

X 

JACOB 

COM2 

XB 

Body  mesh  points 

GRID 

XEX(80, 80 ,1) 

i  =  l,  Cx,  i  =  2,  nx 

JACOB 

COM2 

XEY(80 ,80, I) 

I  =  1,  Cy,  I  =  2,  0y 

JACOB 

COM2 

XMACH 

Moo 

INITIA 

COMl 

XS 

Shock  mesh  points 

GRID 

Y(80, 80) 

y 

JACOB 

COM2 

YB 

Body  mesh  point 

GRID 

YS 

Shock  mesh  point 

GRID 
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3 


A  LISTING  OF  COMPUTER  PROGRAM 


i 

jt 

i 


i 


i 


i 


5 


10 


is 


20 


25 


30 


35 


40 


45 


50 


PROGRAM  NOST IP ( INPUT  *  OUTPUT  ,TAPFS« I NP  JT . T APE 6«0tlTPUT ♦ T APF 1 ♦ T APE  2. 
1TAPE7.TAPE81 

C0MM0N/C0M1/JMAX.KMAX, JM*KM*XMACH*GAM*6AMM) , CN,DT * 5MU. JCS . PpT * 

1  IPRT,H,OM£GA.IT*TAU.ITER.kNT,PTrtRT.PlNF.RINF.OlNF,ClMF*PT.lTS. 

2  IR1 » IN2* 1 AFHD* IGtOM*  TM.IVIS.ITRAN.CF *CC» JNM.RF Y.PRD.CVIS.CV 1  SI • 

3  TWA.ITwA.LIP.KRESfSMJIMP.HTINE .FTINF  »SINF  »EI INF.Rfc YIN*SMM(401 * 
4UETT  (401  *DETL  (40)  *fc  T  (401  .  Th  (401 . 1 TF  *F4CTH,F  ACTT.REYNLfl.PRTURR 

C0MM0N/C0M2/X (40*401 *Y (40*401 *XFX (40* 40*21 • XF Y (40*40*?) ,(1(40.40) 
COMMON/COM3/Q(40*40.41 *EF (40.41 *S (40*40*41 ,5 (4) .At)  (4*4) ,HVFC (40.41 
COMMON/ V I SK/CMUK AP (40 1  * TURMtl (40*401 
C  INITIALIZE  FL0WF1EL0 

CALL  INITIA 

C  PRINT  OUT  STARTING  SOLUTION 

CALL  OUTPUT ( 1 1 
C  OtTFRMlNE  STEP  SIZE 

CALL  EIGEN 

C  COMPUTE  RESIDUE  INFORMATION  AT  START  of  EXECUTION 

CALL  RHS 
CALL  RESIDU 

C  INTEGRATE  EQUATIONS 

DO  1  I*ITS* ITFR 
I T=I T* 1 
CALL  SHOCK 
CALL  0NDRY 
CALL  INTEGH 
CALL  JACOB 

IF (MOO (I.2SI.EQ.01  CALL  EIGEN 
1  CONTINUE 

C  PRINT  OUT  SOLUTION 

CALL  OUTPUT (11 
IK (IV1S.FQ.0)  GO  TO  4 

C  PRINT  OUT  h£AT  TRANSFER  INFORMATION.  ITE*0  WITHOUT  HEAT  TRANSFER 

C  CAL.*  *1  FOR  STANTON  NO.  ONLY.  •?  TEMPERATURE  FIELD  ALSO 

IFdTF.EO.Ol'  GO  TO  4 
CALL  OUTPUT (61 
4  CONTINUE 

C  STORE  STARTING  SOLUTION  FOR  AFTFPBODY  CALCULATION 

C  1 AFBD*1  FOR  STORAGE  OF  STARTING  DATA.  *0  OTrtFRWlSE 

IF ( I AFHO.EQ.O)  GO  TO  3 
CALL  OUTPUT (41 
3  CONTINUE 

C  STORE  SOLUTION  ON  TAPE  FOR  RES1 APT 

C  I w?r 1  FOR  STORAGE  OF  SOLUTION  ON  TAPE?  FOR  RESTART,  *0  OTHEWwUSF 

IF ( IW2.EQ.0)  GO  TO  2 
CALL  OUTPUT (31 
?  CONTINUE 

C  OUTPUT  DETAILED  PESIOJE  INFORMATION 

C  PRINTOUT  LOOP  FOR  K«1 ,KMAX,KRES  AND  FOR  ALL  J.KRES  IS  AN  INPUT 
CALL  RHS 
CALL  RESIDU 
CALL  OUTPUT (51 
STOP 
END 
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SUBROUTINE  ABMATX { J.K ,  I  ) 

COMMON/COM1/JMAX.KMAX,  JM . KM, XMAOH » GAM ,GAMM 1  .CN.OT.SMU.  ICS.PkT. 

1  IPRT.H.OMFGA.IT.TAU.ITFR.FNT.PTORT.PINF .PI  NF ,OlNF .C IMF ,RT « I TS» 

2  IR1.IW2.I  AFBO .  1GEOM.  TM.IVIS.ITPAN.CF.CC.JNM.RFY.PRD.CVI^.CVISI, 

3  Twa.  ITwA.LIP.KRES.SMJIMP.hTINF  .FTiNt  .SlNF,EIINF.RkYlN.Sll«(40)  , 
4L)tTT (40) .OETLI40) .ETI40) ,Th(40) , ! TF ,F ACTh .F ACTT .REYNLD.PBTUPB 

COMM0N/C0M2/X (40.40) .Y (40.40) «XFX (40,40.2) .XF V (40. 40.2) .0(40.40) 
COMMON /COM 3/0 (40.40.4) ,EF (40.4) ,C (40.40 .4) ,& (4) ♦ Ad (4.4) .HVFC  (40.4) 
C0MM0N/C0M4/A (40.4.4) .M40.4.4) ,0(40.4.4) .H0(40.4.4) • 

100(40.4.4) .AX (40) .AX (40) .MX (40) ,PY (40) 

FORM  JACOB  I AN  MATRICES  AT  A  GIVEN  J-h  NODE  POINT.  A  MATRIX  IF  1  =  1. 

B  MATRIX  IF  1=2. 

XX  =  0. 

YY=x£X ( J.K,  I ) 

/Z=XEY ( J.K. I ) 

BI= 1 . 0/Q ( J.K. I ) 

U=0(J.K,2)*RI 
V=G ( J.K ,3) *R I 
SS=GAMM1 *0.5*(U*U»V*V) 

T=YY*U*77*V 

w=GAM*Q(J.K.4)*RI 

AH (1,1) =XX 

AW ( 1 , 2) =YY 

At)  ( 1 ,3)  =  ZZ 

Ab ( 1 .4) =0.0 

Ab (2, 1) =YY»SS-U*T 

AB(2,2)=XX-YY*  < GAM-2.0) *U»T 

AH (?,3) =-YY*GAMMl*V*ZZ*0 

AB ( 2 .4 ) =Y Y*GAMM 1 

Ab ( 3. 1 ) =ZZ*SS-V*T 

AH(3«2)=YY*V-7Z*GAMM1»U 

AB(3,3)=XX-ZZ* (GAM-2. 0 ) *V*T 

Ab (3.4) =ZZ*GAMM1 

AB(4»1)=T*(2.0*SS-W) 

AH (4,2) = (W-SS) *YY-GAMM1*T*U 
Ab (4,3) = (W-SS) *ZZ-GAMM1 *T*V 
Ab (4.4) =XX.GAM*T 
ADO  SOURCE  TERM  IMPLICITLY 

IF(JCS.EO.O.OR.I.EQ.I)  RETURN 
Yl=OT/Y (J.K) 

UD (K , 1 • 1 ) =0.0 
OD(k.1,2)=0.0 
UO  <K ♦ 1 . 3) =Yl 
UD (K  « 1 ,4) =0  •  0 
U0(K.2,1)=-U*V*YI 
UD (K.2,2) =V*YI 
UD (K.2.3) *U*Y I 
OD (K.2.4) =0.0 
tlO  <K,3»  1 )  =~V*V*Y  I 
UO(K, 3. 21=0.0 
UU(K,3»3)=2.0*V*YI 
OD (K ,3,4) =0.0 
UD(K,4.1)=V*(2.0*SS-W)*YI 
UD  <K ,4,2) =-U*V*GAMMl* Y I 
U0(K.4.3)=(W-SS-GAMM1»V*V)*YI 
UD (K ,4.4) =V*GAM*Y I 
RETURN 
END 
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SUBROUTINE  HNDRY 

COMMON/COM1/JMAX,KMAX, JM,KM,XMAru»GAM,(.AMM) ,CN.DT * SMU. JCS.PmT » 

1  lPPT.H,OMEGA,IT.TAU»ITER»fcNT»PTORT.PlNF.RINF.OlNF,ClNF»RT.nS, 

2  Ih1.1W2,IAFBD.IGEOM,TM,1V1S.I TP»N,Cf  ,CC. JN*«RF Y.PRD.CVl S.CVI SI , 

3  TwA,ITwA,LlP*KRfS.SMJIMP,HTINf .FTlNF.STNF.FIINF»RtYlN,SHM(40)  * 
4uETT  (AO)  «DE  TL (40 )  ,L  T  (40)  ,  TH  (40 )  ,  1  Tf  *F  ACTh.F  AC.TT  »R£YNLD,PRTltRR 

C0MM0N/C0M2/X (40*40 ) ,Y (40*40) ,XFX (40.40.<?> «XF Y(40,4(>»?> ,0(40,40) 
COMMON/C 0M3/0 (40,40,4) ,EF (40 ,4) ,C (40,40 ,4 ) « G (4) , AS (4,4) .HVFC (40,4) 
DIMENSION  P (40 ,3) ,PXI (40) .PET A (40) ,U{40,3> ,UXI (40) .UETAUO)  . 

1 V (40 , 3) , VX I (40) ,  VE TA  (40 1 , P (40 , J ) 

DIMENSION  T (40,3) ,bOIAG(40) ,01 AO (40), AO I AG (40) ,PIGHT(40) 

DIMENSION  DUMMY (40) ,DCON(40) ,tCOM(40) 

...  DATA  Cl, C2.C3/-3. 0,4.0, -1.0/  THIS  SFT  DATA  USED  FOP  3  POINT 

...  ONE  SIOE3  DERIVATIVE  APPROXIMATION  AT  BOOY  BOUNDARY 

...  DATA  Cl, C2,C3/-2. 0,2.0, -0.0/  THIS  SET  DATA  USED  FOw  ?  POINT 

...  ONE  SIDED  DERIVATIVE  APPROXIMATION  AT  BODY  BOUNDARY 

DATA  Cl, C2.C3/-3. 0,4.0, -1.0/ 

...USE  REFLECTION  TO  SIMULATE  PLANE  OE  SYMWFTRY  AT  J=2 
DO  12  K= 1 »KMAX 
(J(1,K,1)=Q(2,K,I) 

0(1,K,2)=Q(2,K,2) 

0  ( 1  ,K ,4) =Q (2,K,4) 

12  Q(l»K«3)=-Q(2,K«3) 

...USE  FIRST  ORDER  EXPTAPOLATION  TO  SIMULATE  SUPEPSONIC  OUTFLOW 
...BOUNDARY  CONDITION  AT  JMAX 
DO  1  N-l ,4 
DO  1  K=  1  , KM 

1  0(JMAX,K,N)=(2.0*0( JM,K,N)-0( JM-1 ,K,N) ) 

IF(JVIS.EQ.O)GOT014 

...APPLY  VISCOUS  NOSLIP  BOUNDARY  CONDITION 

ITWA=0  FOR  ADIABATIC  WALL  *  ITWAsl  FOR  ISOTHFPMAL  WALL 
DO  15  J=1,JMAX 

DCON(J)=XEX(J,1,1)*XEX(J,1,2)*XFY(J,1,1)*xEV(J,I,2) 

ECON(J)=XEX(J,l,2)*«2*XEY(J*l»2)**2 

DO  15  K=l,3 

P(J,K)=(0(J.K,4)-0.5*(0(J»K,2)**?«QCJ,K,3)*»2)/0( J,K,1) ) 

>  *0(J.K)*GAMMl 
T(J,K)=P(J.K)/0(J,K)/3(J,K,l) 

15  CONTINUE. 

...SET  UP  COEFFICIENT  MATRIX  FOR  TRIOTAGONAl  INVERSION 
00  16  J=?,JM 
bDI AG( J) =-OCON ( J) 

UIAG(J)=CI*ECON( J) 

ADI  AG (J) sOCON ( J) 

16  CONTINUE 

0 1  AG (2) =01 AG (2) -OCON (2) 

...COMPUTE  WALL  PRESSURE 
DO  17  J=2,JM 

17  RIGHT ( J) =EC0N( J) * (-C3*P ( J, 3) -C2*P ( J,2) ) 

RIGHT ( JM) =RIGHT ( JM) -OCON ( JM) *P (UMAX , 1 ) 

CALL  TRI8(BDIAG,0IAG,ADIAG,DUMMY,RIGHT,2,JM) 

DO  18  J=2.JM 

18  P ( J, 1 ) =RIGHT ( J) 

P(l,l)-P(2,l> 

...COMPUTE  WALL  TEMPERATURE  FOR  ADIHATIC  WALL 
If  (ITWA.FQ.l)  GO  TO  21 
00  IP  J=2.JM 

19  RlGHT(J)sECON(J)*(-C3»T(J,3)-C2*T(J,2> ) 

RIGHT (JM)=HIGMT ( JM) -OCONI JM) »T (JMAX. 1 ) 

CALL  TRIB (HOI  AG, 01  AG, AO  I  AG, DUMMY, RIGHT, 2, JM) 

DO  20  J=2,JM 

20  T(J,1)=RIGHT(J) 

T(l,l)*T(2,l) 

21  CONTINUE 

...COMPUTE  CONSERVATIVE  VARIABLES 
00  22  J* 1 , JMAX 
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IF  (ITWA.EQ.O)  T*A=T<J,1> 

T 1 =TWA 
Pl=P(J.l) 

70  R1=P1/T1 

01=1. 0/D ( J. 1 ) 

0<J,1,1)=R1*0I 
Q(J, 1(2)30.0 
Q( J, 1.3)30.0 

75  q<j,i,4)=P17GAMM1*DI 

22  CONTINUE 
RETURN 
14  CONTINUE 

C.. .APPLY  INVISCIO  BOUNDARY  CONDITION 
80  C.. .SATISFY  TANGENCY  CONDITION  USING  CHARACTERISTIC  EQUATION 

UO  3  K=1 *3 
DO  3  J*l .  JMAX 
7  =  1 .0/0 ( J(K ( 1  ) 

R<J.K>=0(J.K,1>*0(J.K) 

85  U(J,K)=0(J(Ki?)*2 

V(J,K>=Q<J,K,3)«Z 
t2--0(JiKi4l‘0IJ.Ki 

3  P<J.K)=(E2-0.5*R< J»K)*(U(J.K)**?.V( J. <)**?) )*GAMMl 

C.. .COMPUTE  P-XI,  U-X I • V-X I (  P-ETA(U-FTA(  AND  V-ETA  DERIVATIVES 
90  DO  4  J=2»JM 

PXI (J)3<P(j*l,l)-P(J-l,l) )«0.5 
UXKJ)  =  IUIJ*1.1)-U(J-1(1)  >*0.5 

4  VXI <J)=<V<J*1.1)-V<J-1,1> )*0.5 
PXI < 1 ) =-PX I (2) 

95  UXI (1)=-UXI (2) 

VXI <1>=VXI (2) 

PXI <JMAX)  =  (3.0*P( JMAX, 1)-4.0*P(  |M.1>.P<JM-1.1>  >*0.5 
UXI ( JMAX)=(3.0*U<JMAX.1)-4.0*U(JM(1 >.J(JM-1.1 ) )*O.S 
VXI <JMAX)=<3.0*V<JMAX,l)-4.0*V<JM.l).v< JM-1.1) >*0.5 
100  DO  5  J= 1 , JMAX 

PETA( J>=<-3.0*P<J»1>*4.0*P<J(2)-P<J.3> )*0.S 
uETAt J)=t-3.0*U(J*l>*4.0*U(J.2)-|l(J.3> )*0.5 
VtTA< J)s(-3.0*V<J.l)*4.0*V(J,2>-V(J.3> )*0.5 

5  CONTINUE 

105  K  =  1 

IF(IT.FQ.ITER)  wPI TE ( 6. 1 02) 

102  FOHMA75 *OFROM  SUB.  BNDRY* ) 

DO  2  J= 1 ( JMAX 

CBb=SOPT(GAM*P» J.l )/R( J.l ) ) 

110  7=1 . 0/SQRT  <XEX(J.1»2) **2*XEY ( J. I ,2) **2> 

UBARslM  J(  1 )  *XEX  ( J«  1 ♦ 1 ) »V< J.  1 )*XEY< J(  Id) 

VBAR=U (J,1)*XEX(J»1»2)*V(J,1)*XEY(J.1»?) 

F  F.=U8AR*PX  I  (J)  *R(  J.1>*CHR**2*<AFY<J»  1,1)  *11X1  < J) »XE Y C J, ) . 1 ) *VX I  (J)  ) 

>  -CflH'Z* (XtX  < J, 1,2) *< XIX (J.l ,1 ) *PXI ( J) ♦&< J, 1 ) •U8AR*UXI ( J) ) ♦ 

115  >  XEY(J(1.2)*(XEy(J<1.1) *PX I(J)+p(J,1) *UrtAp*VX I ( J) ) ) 

1  ♦CVIJ(1)/V(J.D) 

HTAU=CPH/Z*PETA(J>-R{J,1)*CRR**?*(XEX(J,1.2>*UETA(J>*XFY<J.1.2>* 

>  VE  TA  I  J)  )  -EE 
P1=P<J.1).PTAU*OT*0.2 

120  IF(Pl.LF.O.O)  GO  TO  9 

10  CONTINUE 

IF ( J.GT.2)  GO  TO  11 

A=(P(4»l)-25.0*P(3»l)/9.0*16.0*PT/9.0)/6.?5 
H=(P(3»l)-PT-27,0*A/8.0)*4.0/9.ft 
125  Pl=PTtA/8.0*d/4.0 

11  CONTINUE 

Rl=(Pl/ENT>**( 1.0/GAM) 

01 =SQRT (2.0*GAM/GAmmI *ABS (PTORT-Pl/Rl ) ) 

IF ( ABS (XEY (J.l. 2) ) -0.000001)  6.6,7 
130  6  TMtTh=l. 570796327 

GO  TO  8 

7  THF. TB=ATAN(-XEX(J»1»2)/XEY(J»1»?)  ) 

8  CONTINUE 

Ul  =01 *COS ( THETB ) 

135  Vl=01*SlN(THETB) 
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TH£TBD=90.0-TH£TB«S7.2967H' 

IF(IT.EO.ITER) 

>WRIT£(6,101)  J.C8B»Z»PTAU.DT.XtX(J,l,l),xFY(J.l.l>.XEX(J,1.2)  . 

>  XEy( J,1 ,2) ,PXI < J) .PETA ( J) »P( J.l ) .PI ,UXI ( J) ,IIFTA( J) ,R< J, ) ) ,R1 , 

>  VXI (J) , VETA ( J) ,U(J.1),U1»V(J,1>,V1»THETH0»EF 

101  FORMAT!*  J=*,I2,4X.BE13.5,/,9X,PF)3.5,/.9X,PF13.5> 
DI=l.fl/0(J,l) 

(J(J,1.1)=R1*DI 

0<J.1.2)*R1*U1*OI 

Q(J»1«3)=R1*V1*DI 

Q< J.l »4>s(P1/GAMM1*0.5*R1*Q1**2>*OI 
2  CONTINUE 

he Turn 

9  WRI T£ (6. 100)  IT.J.Pl 
P1=ABS(P1 ) 

GO  TO  10 

100  FORMAT!*  IT£R=*,I4.3X,*J=*,I2,3x,*P1**,F12.4) 

END 


SUBROUTINE  BTRI(IL.IU) 

C0MM0N/C0M2/X (40.40) , Y (40 .40 ) . XF*  <40 » 40 .2) , XF Y (40 .40.2) .n (40 .40) 
C0MM0N/C0M3/Q  (40.40.4)  «EF(40.4)  .S<40.40*4)  .0(4)  .AK4.4)  .HVFC(40.4) 
COMMON/COM4/A (40.4.4) ,b(40»4.4) .0(40.4.4) «h0(4O.4.4) . 

1(10(40.4.4) .AX (40) .AY (40) »HX(40) ,PY(40) 

COMMON/LU0/  Lll  .L21.L22.L31.L32.L33.L4) .L42.I  43.L44.V1 , V2.V3.V4, 

1  U12.U13.U14.U23.U24.U34 

DIMENSION  H(4,4) ,F (40,4) 

REAL  L11,L21.L22.L31.L32.L33.L41,L42,L43,L44 

EQUIVALENCE  (EF ( 1  .  1 ) , F  ( 1  ,  1 )  ) 

C  INVFRSION  OF  BLOCK  THI 01  AGONAL .. .  A.B.C  ARE  4*4  Rl OCKS 
C  F  IS  FORCING  FUNCTION  AND  SOLUTION  IS  OUTPUT  IN  F.  3  IS  OVERLOADED 
C  BLOCK  INVERSIONS  USE  N0N2I VOTED  LU  DECOMPOSITION 
C  IL  AND  IU  ARE  STARTING  AND  FINISHING  INDICES 
IS  =  IL  *1 
I  =  IL 
DO  11  N=l,4 
DO  11  M=l,4 

11  H(n,m)=B(I.N,M) 

CALL  LUDEC(H) 

01  =  Vl*F (  I  .  1  ) 

02  =  V2* (  F ( I  ,2)  -  L21 *01 ) 

03  =  V3*(  F ( I ,3)  -  L31*01  -  L32*D2) 

04  s  V4*<  F  ( 1 ,4)  -  L41*fJl  -  L42*D2  -  L43*03) 

F (1,4)  s  04 

F (1,3)  =  03  -  U34*D4 

F ( 1 » 2 1  =  02  -  U24*D4  -  U23«F (1.3) 

F ( I , 1 )  =  01  -  U14*D4  -  Ul 3*F (1.3)  -  U12*F(I,2) 

DO  12  M= 1,4 

01  =  V1*C ( I . 1 *M) 

D2  =  V2*(  C ( I »2»M)  -  L21 *01 ) 

03  =  V3*(  CII.3.M)  -  u3 1 *01  -  L3?*D2) 

04  =  V4*(  C(I.4,M)  -L41 *01  -  L4?*D2  -  L43*03) 

B ( I .4  »M)  =  04 

B < I , 3.M)  =  03  -  U34*04 

H  ( I » 2, M)  =  02  -  U24*D4  -  |t?3*B  ( 1 ,3.H) 

B(I.l.M)  =  01  -  Ul4*04  -  U13*B ( I ,3, Ml  -  U12*6 ( I »2»M) 

c. . .forvaro  sweep 

12  CONTINUE 

no  i3  isis.iu 
1R  =  I  -1 
DO  14  N=1.4 

F  ( I «N)  =  F(I.N)  -A  < I , N. I ) *F ( I«* 1 )  -A ( I ,N,2) *F ( IP, 2) 

1  -A(I,N,3)*F(IR.3)  -  All, N,4)*F (19.4) 

00  14  M= 1 *4 

H (N.M)  =  B(I.N.M)  -A(I,N,1)*B(IR.1.M)  -A(I.N.2)*B(IR.2,M) 

1  -A(I,N.3)*B(IR.3,M)  -  A(I»N,4)*B(IR.4,m) 

14  CONTINUE 
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60 
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75 


80 


CALL 
01  = 

02  = 

03  = 

04  = 

F (1,4) 

F (1,3) 

F ( I ,2) 

F  ( I  ,  1 ) 
IF  (  I  - 


LUOEC(H) 

V1*F  < I »  1 ) 

V2*(  F ( I ,2) 

V3*(  F ( I ,3) 

V4*(  F ( I ,4 ) 

=  04 

=  03  -  U34*04 
=  02  -  U2A*04 
=  01  -  U 14*04 
IU) 16,13,13 


-  L21*D1) 

-  L31*01  -  L 32*02) 

-  L41*01  -  L42*02  -  L43*03) 


U?3*F (1,3) 

01 3*F (1,3)  - 


U1?*F (1,2) 


16  CONTINUE 
00  15  M= 1 ,4 
01  =  V1*C(I,1,M) 

02  =  V2*<  C(I,2,M) 
03  =  V3*<  C(I,3*M) 
04  =  V4*(  C(I,4,M) 
Ft  ( I  ,4»M)  =  04 
ti(I,3,M)  =  03 


-  L21*0l) 

-  L31*01  -  L3?*02) 
-L 4 1*01  -  L4?*l)2  - 


L43*D3) 


U34*D4 

h(I,2,M)  =  02  -  U24*D4  -  U23*B ( I »3*M) 

0(1,1 ,M)  =  01  -  U14*D4  -  U13*B(f,3,M)  -  012*B(I,2,M) 

15  CONTINUE 
13  CONTINUE 

C...bALK  substitution 

IT  =  IL  ♦  IU 
00  21  II  =  IS,IU 
1  =  IT  -  II 
IF*  =  1*1 
00  22  N=l,4 

F(I.N)  =  F  < I *N)  -  B ( I  ,N, 1 ) *F ( IP, 1 )  -  B(I  ,N,2) *F ( IP, 2) 
1  -B  < I  ,N,3)*F (IP, 3)  -  B ( I  ,N,4)*F (IP, 4) 

2?  CONTINUE 
21  CONTINUE 
WE TURN 
FNU 


.^1 


I 


10 


20 


?5 


30 
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SUBROUTINE  CTXY(XAA,YAA,DtD.HF,CT> 

C0MM0N/H0TH/X01  , X02  ,  X  03  »  X04  ,  Y01  ,  Y02  ,  V03  ,  T04  ,  Si  »  .si 
1W3.W4.CT1 ,CT2,CT3,CT4,CT5,CT6,xnn,Wd03Y 

C0MM0N/XYPS/X1 ,x2*X3,X4,X5,X6,X7,Y1,Y2,Y3,Y4.vc.v*.v- 

THIS  SUHROUTINE  CALCULATES  X»Y,WF  AND  OFLTA  f^d  a  r*VFi 
P2=2 •* AT  AN( 1 , ) 

IF (CT.GT.CT1)  GO  TO  2 
A=XOO*TAN(CT) 

0=1 »  *TAN (CT ) **2 
H=(X01*A*TAN(CT) )/0 
C= (X0l**2*A**2-Rl**2> /O 
X=b-SORT (b**2-C) 

IF (Rl.GT.10.)  X=0. 

GO  TO  10 

IF (CT.GT.CT2)  GO  TO  3 
A=X00*TAN(CT)-Y0? 

0= 1 , ♦ TAN (CT ) **2 
H=(X02*A*TAN(CT) >/Q 
C=(X02**2*A**2-R2**2)/Q 
*=b-SQHT  <B*»2-C) 

’0  10 

IF (CT.GT.CT3)  GO  TO  4 
A1  =  XOO*T  AN (CT ) -Y2 

X=(X?*A1*TAN(SL1 ) )/(l.*TAN(CT)*TAN(SLl) ) 

<>0  TO  10 

IF (CT.GT.CT4)  GO  TO  5 

A=XOO*TAN (CT ) -Y03 

0=1.«TAN(CT)**2 

H= (X03*A*TAN(CT) )/0 

C=(X03**2*A**?-H3*»2)/0 

X=B*SQRT(B**2-C) 

GO  TO  10 
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5  If (CT.GT.CT5)  GO  TO  6 
A1=X00*TAN(CT)-Y4 

X=(X4.A1*TAN(SL2>  )/(l.*TAN(CT)*TAN(SL2)  > 
GO  TO  10 

6  1F(CT.EG.P2>  GO  TO  11 
A=X00*TAN(CT)-Y04 
0=1.*TAN<CT)**2 
Ei=(X04*A*TAN(CT)  )/Q 

C=  (X04**2.A*«2-R4**2> /Q 
X=B-SQRT(H**2-C) 

10  HF  = ( XOO-X ) /COS  <CT ) 

OE  D=XOO-BF 
Y=(XOO-X)*TAN(CT> 

GO  TO  12 

11  X=XOO 
Y=RBODY 
DEl)=XOO-Y 
hF  =  Y 

12  CONTINUE 
X  AA=X 
YAA  =  Y 

Rt  turn 
end 


SUBROUTINE  OISSIP 

COMMON/COM 1/JMAX.KMAX, JM.KM.XMACH.GAM.GAMM1 ,(  N,OT»SMU,.lCS.PBT. 


1  IPRT.H.OMEGA,  IT.TAU.IIFR.ENT.E'TOkT.PIMF.BInF ,0  INF ,C INE .PT, I TS, 

2  IHl.IW2.IAFBn,IGE0M.TM,IVIS» lTPXN.CE ,CC, JNM.PFY.RRD.CVlS.CVlSl , 

3  TWA.ITWA.LIP.KRES.SMJIMP.HTINE  .FTINF.MNF.EIINF.REYIN.SIMUO)  . 
40ETTI40) .OETL (40) .E T ( 40 ) .TH<40 ) , T TF »F AC TH ,F ACTT .RE YNLD.PPTIIPB 

COMMON/COM2/X (40.40) • Y (40 .40 > . XFX (40 . 40 • 2 > . XE Y (40 .40 . 2 > • 0 (40 . 40 ) 
C0MM0N/C0M3/Q (40.40.4) .EF (40 .4 ) ,S (40 . 40 .4 ) , G (4) • AB ( 4.4 ) »HVF C (40,4 ) 
SMOOTH  IN  THE  X  ANO  Y  DIRECTIONS  AMU  ADD  SMOOTHING  TtRM  TO  S  ARRAY 
OATA  C1S»C2S,C3S,C4S/-2.0»S.0, -4.0, 1.0/FOR  LTNFAR  t X TRAP  AT  SHOCK 
DATA  ClS»C2S»C3S.C4S/-1.0»3.0»-3.U»l,0/EOR  PaRAB  ExTPAp  AT  SHOCK 
OATA  ClB,C2B»C38»C4B/-2.0»S.0»-4.0»1.0/EOR  LINEAR  txTPAP  AT  BODY 
DATA  C18.C 2B, C 3B.C4B/-1. 0.3. 0.-3, 0.1.0/F UR  PARAH  tXTRAP  AT  HOOY 
DATA  C 10,C20,C30»C40/-2.0,S.O ,”4 .0,1.0/FUP  LINEAR  EXTPAP  AT  OUTFI.O 
DATA  CIO, C20.C30.C40/-1. 0,3. 0,-3. 0,1. 0/FoR  PAPAB  EXTRAP  AT  OUTElOW 
DATA  C1S.C2S.C3S.C4S/-1. 0.3. 0,-3. 0,1.0/ 

DATA  C1B.C2B.C38.C4B/- 1.0. 3. 0,-3. 0,1.0/ 

DATA  C10.C20.C30.C40/-2. 0,5. 0,-4. 0,1.0/ 


KMMsKM- 1 
JMM=JM-1 

..SMOOTHING  IN  XI  DIRECTION 
DO  4  N=1 ,4 
DO  2  K=?,KM 

..USE  LINEAR  OR  PARABOLIC  EXTRAPOLATION  FOP  J=JM 
..SEE  DATA  STATEMENTS  ABOVE 

S(JM,K»N)=S(JM,K»N)-SMU*0.125«(C10*Q(JMAX,K.N)*D(JMAX,K) » 

>  C20*Q< JM,K»N)*D< JM.K) .C30*Q< JMM.K.N) *D  < JMM.K) ♦ 

»  C40*Q(JM-2»K.N)*D( JM-2.K) )/D(JM,K> 

DO  2  J=3,JMM 

2  S(J.K»N)=S( J.K,N)-SMU*0. 12B* (Q ( J-2.K » N> *D ( J-2.K > -4. 0*0 ( J-l .K.N> 


>  *0(J-1,K) .6,0*0 ( J.K»N)*0 ( J.K) -4.0*0 ( J.l.K.N) *D(J*1 ,K) *0( J« 2, K,.: 

>  *D(J.2,K) )/D(J,K) 

..SMOOTHING  IN  ETA  DIRECTION 

DO  1  J*2, JM 

..USE  LINEAR  OR  PARABOLIC  EXTRAPOLATION  AT  BODY  AND  SHOCK 
..SEE  OATA  STATEMFNTS  ABOVE 

S(J.2,N)=S<J»2»N)-SMU*0.125*(C1B*0(J.1,N)*0(J.I). 

>  C2B*Q(J,2,N)*D(J,2) *C3B*0 ( J.3.N) *0 ( J» 3) ♦ 

>  C4B*0 (J,4»N)*D(J»4) )/D(J,2) 

S( J,KM,N)*S( J.KM»N)-SMU*0. 12S* (f I S*Q( J,KMAX.N)*D( J.KMAX) ♦ 

>  C2S*0( J,KM.N)*0( J.KM) .C3S*0( J.KMM.N) *0 ( J.KMM) ♦ 

>  C4S*Q(J.KM-2,N)*D( J.KM-2) )/0< J.KM) 
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00  1  K=3,KMM 

1  S(J,K.N)=S(J.K,N>-SMU»0.125*(Q(  I, K-2.  St )  *0  ( J, K-? ) -4 . 0*0  ( J.K- 1  , N) 

*5  >  *0<J.K-1).6.0*Q<J.K.V)*D<J,K>-4.0«0(J,K.1,N)*0(J.K.1).0(J,K.2.N> 

>  »U|J,K.2)I/0(J,K) 

4  CONTINUE 

C.. .COMPUTE  SMOOTHING  FOR  J  =  ?  RY  USING  SYMMETRY  CONDITIONS 
1)0  3  K  =  2  *  KM 

5 0  S(2,K.1)=S(2.K,1> — SMU*0 • 1 25*  ( “4 .0*0(1. K.1)*D(1«K) *6.0*0 (2.K.1)* 

>  D(2,K)-3.0*Q(3,K.l>*3(3.K).a<4,K.l)*3(4.K) )/0(?.K) 
S(2,K.2)=S(2.K.2)-SMU*0.125*(-4.n*Q(l,K.2)*:)(l,K)*6.0*0(?.K,2>* 

>  L>(2.K)-3.0*Q<3.K.2)*3(3.K).O<4,k,2)*D(4.K))/D(?.K> 

S  (2 ,K *3)  =S ( ?,K » 3)  -SMU*0. 1 25*  (-*». p*o  ( 1  .k  * 3)  *D  ( 1  .K )  ♦{>•  0*0  ( ? *K *  3)  * 
55  >  D<?,K>-5.0*u(3.K.3)*3(3.K).0(4,k.3)*3(4,k))/1)(?.<) 

S<2,K.4)=S(2.K,4> -SMU*0. 125* (-4.0*U( 1 .* .4) *3 ( 1 ,K ) *6. 0*0 ( ?«K ,4 ) * 

>  L><2.K)-3.0*G<3.K.4)*3(3.K) .o(4,K.4)*D(4.K) )/D(2,K) 

3  CONTINUE 

ML  T  URN 

60  L  NO 
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subroutine  feconij.k.I) 

CUMMON/COMl/JMAX ,KMAX, JM,KM.XMArH,GAM,&AMMl .CN.DT  »5MU.  ICS.  PR  T ♦ 

1  IPRT.H.OMEGA.IT.TAU.ITFR.f NT.PTORT.PINF.Pplf .OINF.CINF.PT.ITS, 

2  Ik! . IW2. I  AErifl, IOtOM, TM. I VIS, I  I  DA N , CF . CC . JN« . PE Y . PRO . C V I  5 . C V I S 1 . 

3  TwA.ITwA.L  IP.KRf  S.SMJIMP.HT  INE  ,  F  T  I  NF  ,  S  I NF  ,  E  I  I  E'F  »  Rt  Y 1 N  ,  Sum  ( 40  >  . 
4DETT (40) ,OFTL (40) #t  T (40) .TH(40) . TTE ,F ACTri  ,F ACTT .RE YNLD. PPTURB 

COMM0N/C0M2/X (40*40) »  Y (40*40) *  XF  X (40.40.2) »  XF  Y (40*40,2) ,0 (40.40) 
CUMM0N/COM3/0 (40*40,4) ,EF (40  * 4) . « (40  *  40 .4 ) * G ( 4 1 . AB ( 4 .4 ) .HVEC (40,4) 
DATA  MVFC/160«0.0/ 

C . .  ,FOljM  e  CONSERVATIVE  VAR  I  AHLtS  (1  =  1)  OR  F  CONSERVATIVE  VAPIABl.ES 
C. , . (1=2)  AT  A  GIVEN  NOOt  POINT 
w=Q(J.K, 1 ) 

Pl=l .0/w 

U=0(J*K,?)*RI 

V=Ci(J»K,3)*RI 

P0J=GAMM1*(Q( J.K.4 >-W*0.5* (U*U*V*V) ) 

XX  =  0. 

YY=XFX ( J.K. I ) 

Z=XFY(J,K,1> 

CAPUV=XX*YY*U*Z*V 
G( 1) =W*CAPUV 

G(2)=Q< J*K.2)*CAPUV*YY*P0J 
G ( J ) =Q ( J ,K , 3) *CAPUV*Z*P0J 
G  (4)  =0(  J.K.4)  *CAPUV*  (CAPUV-XX)  «PGJ 
C.. .SOURCE  TERM  IN  ETA-MOM.  tON.  FOR  AXlSYMMETRIC  FLOW 
IE (JCS.EO.O.OR.I.EQ.l)RETURN 
Y I  =DT/Y (J.K) 

HVtC(K.l)=Q(J,K,3)*Yl 

HVtC (K*2) =0( J.K, 3) *YI*U 

HVtC  (K.3)  =()(  J.K. 3 1  *YI*V 

HVEC (K.4) = (0( J.K.4) *POJ) *V* Y I 

RETURN 

ENU 


1  SUBROUTINE  EIGEN 

CUMMON/COMI/JMAX.KMAX,  JM.KM.XMACH.GAM.GAMM)  ,cn.(">t«smu..ic<?,prt. 

1  1PPT  .H.OMFGA,  I  T.IAU,  ITFR.ENT.PTDRT.PIE'F.PINF.OIN'  .ClNF.l  T.  ITS, 

2  IR| . IW2, IAFHD, IGEOm, TM, I  VIS. I IRAN.CF ,CC . JNM, wF Y.PRD.CVI « ,C V I S I , 

5  3  TwA.ITWA.LlP.KHf 5, SM JI MR ,HT I NE  , F  T I NE  .SINE  .f I  1 NF • Rt Y I N. St  >M 1 40 ) . 

4UETT (40) .Of  TL(40i * E T (40) »TH(40) .TTE .E AT Th , F  AC T T . RE YNL D , PPTURH 
C0MM0N/C0M2/X (40.40)  . Y (40,40 ) ,XE  X (40.40.^1 .XE  Y( 40. 40, 2) .0(40.40) 
C0MM0N/C0M3/N (40 « 40 . 4 ) ,EF  <40 , 4 ) , S (40 • 40 .4 ) , 3 ( 4 ) . A3 ( 4 .4 ) .HVE  C (40.4 ) 
C.. .COMPUTE  STEPSIZE  GIVEN  COURANT  NumrE r 
10  IF ( 1PRT.GT.0)  WRITE  (G,:n0) 

S I  <»M  A  X  =  0  •  0 
SIGMIN=10.E*100 
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no  i  k= 1 *kmax 
uo  i  j= l «  jmax 

kI=1.0/0<J,K,1) 

U=Q(J.K,2)*»I 
V  =  (J(J.K,3)*RI 

XX=GAM«GAMMl<*(0<  J.K.4)  «RI-0.5"  (ll*U*V«V)  ) 

IF  <  XX  >  2,2,3 

?  WHITE (fc.103)  J.K.Ol  J.X.4)  .HI .U.V.XX 
XX=-XX 

3  SPSNO=SOHT < XX) 
xix=xfX(j,K.n 

XIY=XEY(J,K.1> 

ETAX=XEX(J.K,2) 

ETAY=X£Y(J,K,2) 

XE  T  =  0  « 

SlGA  =  ABS(XE  T*U*X1X*V»XIY> ♦SPSND«S(jHT  <XIX**2«XIY**?> 

SIGfi=AHS  EXE T*U*ETAX« V*£TAY) ♦SPSRn»3QHT (t 7  Ax »*?»FT4 Y*»2> 

SIGAH=AmAX1  (SIGA.SlGb) 

S I GABM- AM  INI ( S IGA  »  S 1 G3 ) 

IF <SIGAb.GT.SIGMAX)G0T04 
GOTOS 

4  SIgmAX=SIGAH 
Jt IGMX=J 
KtlGMX=K 

5  CONTINUE 

IF (SIGAHM.LT.SIGM  IN  1007  00 
GOTO? 

6  SIGMIN=SIGABM 
JE 1GMN=J 
KF.lr,MN  =  K 

7  CONTINUE 

IF  ( IPHT.EO.O)  GO  TO  1 
WHITE  (6.101)  J  ,  K  »  S  l GA  ,  S I GH 
1  CONTINUE 

(>T=CN/AHS(SI  GMAX ) 

WHITE  (6. 102) SIGMA X, JEIGMX.KF I GMX , S I GM I N *  JF I GMN.KE IGMN.CN.OT 
n  =  O.N<*OT 

100  FORMAT  (<*0'»»3X,»J*.4X,»K*, ?X . *S 1 GA * . BX . *S I GR * ) 

101  FOkmAT (2IS.2F12.6) 

>  U?  F  OrmAT («0»,«>SIGmAX  =  »,E10.4.3X,*  l  =  *  »  I  b  .  3X  ,  *K  =  *  , l 5 . 3A . *S I Gm  I  N  =  * . 

>F  .  . !.3x,*j=«,IS,3X,*K=*,IS,3X,«rN=*,E10.4,3X.*DT=«.£10.4) 

10.1  F  rma,  (*(,Nf  GATIVF  SORT  IN  EIGEN  at  J=  *  .  J  2  *  *  K  =  » , ] ? , 3 X , *E / J=* , E ) n  .  4 
>  tlx . «J/H=*,E10.4,JX,*U=“»F l0.4.1X,»V=®.tl0.4,3X.®UISCRM=*,E 10.4) 
hF Town 
end 


SUHROUTINE  F  TATME  T.C.  ,KMAX  ) 

'IMF  NS  I  i  N  JJI (3) » JJF  ( 3) .XXI  ( 3)  »XXF  ( J ) .00  *  I  ( 3 ) . DOXF ( 3 ) . F T (40 ) 
i  *U  A  JJI  (  1  ) . JJI  (?)  ,  JJI  (  T) .JJF ( l )  , JJF (2) .JJF (3) /I . 1 B ,43. 1 7 ,4? , 4H/ 
l»«  '  X  XX  I  (  1  )  .XX  I  (?) ,XXI (3) /O. ,0. Ilb.0.64/ 

I*  A  T  T  XXF  (  1  )  ,XXf  (?)  .XXE  (3) /(I.  1 ,0.6  .  J  ./  ' 

OA  (  a  u()x  I  ( 1 )  .00X1  (2)  .  30X  I  (  J)  /O.nnl  .O.OIS.O.or/ 
ija r a  unxF ( i ) .nuxF  (?) . joxf <3)/o.nib,o.o3»o.ous/ 

.Ml  =  KMAX-1 

rat  =  (CF» I . ) / (CF-1 . ) 

(il  TAC  =  1  ./KM1 
.7(1)  =  0. 

<•  T  (KMAX  1  =  1. 

DO  1  X  =  2, KM  1 
F  TAC  =  IK-1 > »0E  TAC 
t  X  =  l.-FTAC 
AHO  =  RAT**EX 

I  F  7  ( K )  =  1.  ♦  CF • ( 1 , -AHG) / ( 1 . ♦ AHG) 

HE  TURN 
END 
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SUBROUTINE  GRIO 

COMMON/COM 1 / JMA* * KMAX t JM, KM, XMArHtOAM.OAMMl ,TN * OT * SMU. JC« *PRT • 

1  IPRT  *H*OME  GA  « IT«TAU*ITER»tNT*HTORT*RlNF«PINK*QlNF»CINF*PT*ITS* 

2  IHI » IW2* I AFHO  * IGEOM*  TM*IVIS*IlPAN*CF  *CC» JNM.RF Yt PRD*CvlB*CVI51 • 

3  TWA* ITWA.L IP tKRFS*  SMJIMP.HT INF ,FTINF  *SINF*E1  INF.REYIN.SIIM(40>  . 
4DETT(40) *DETL(40).ET(40) . TH (40) , I TF *F ACTM.FACTT.REYNLD.PPTURB 

COMMON/COM2/X (40*40) . Y (40*40 )« XFX (40.40.2) »XE Y (40*40 ,2) *0(40*40) 
C0HM0N/C0M3/Q (40*40*4) *FF(40*4) *S (40.40 .4 ) *G (4) *Ad(4«4) *HWFC(40,4) 
C0MM0N/BDTH/X01 .X02.X03.X04.Y01 ,Y02. Y03.YO4.SL1 .SL2.SL3.P1 *R2« 
lR3.R4»CTl.CT2.CT3.CT4.CTS.CT6.Xnn.RttG3Y 
C  THIS  SUBROUTINE  Dt TERMINES  X  ANO  Y  FOR  6RI0  POINTS 
Pl2=2.*ATAN( 1 • ) 

PI=2.*PI2 

OTR=PI/lflO. 

TM=(90.-TM)*DTR 

JNM1 = JNM- 1 

TH( JNM) =  TM 

JNMS=JNM*ITRAN 

DTH= (PI2-TM) /FLOAT ( ITRAN) 

JNMP= JNM* 1 
TMM=PI?-TH 

C  FOR  SPHERE  PORTION 

UTHMI N=TM/ (FLOAT (JNMl) -0.5) 

UTH1=0.5*DTHM!N 
00  SI  J=1.JNM1 
51  TH(J)*(J-l)*DTHMlN-DTHl 

IGEOM*0  UNIFOPM  POINTS  ON  SPHEKF  FOR  SPHtRF-CONF 
IGE0M41  READ  IN  XB.YB.XS.YS 

IGE0M=2  READ  IN  TH(J)  ANO  OETT(J)  FOR  ARBITRARY  BODY  SHAPE 
IGE0MS3  UNIFORM  SPACING  FOR  THU)*  CAL.  OFTT(J)  AND  DETFRMInE  THE 
XB  ANO  YB  FOR  INOENTEO  NOSE 

IGE0MS4  READ  IN  TH(J>  ON  NOSE  ANO  CAL.  OtTT(J)  AND  FINAL  XB  ANO  YB 
IF(IGEOM.EQ.O)  GO  TO  41 
IF ( IGEOM.EO. 1 )  GO  TO  3 
IF ( IGEOM-3)  42.43.43 

42  REAO(5*102) (TH(J) *J=1*JMAX) 

RE  AD  (5*  1 02)  (0ETTU>*J*1*JMAX) 

102  FORMAT (8F10.5) 

WRITE (6*103) 

103  FORMAT (/.2X»«REA0  IN  TH(J>  ANO  OFTT (J) *./.?0X ,«J*.5X**TH( J) .DEGREE 
1*.5X.*DETT(J)*> 

DO  104  Jsl.JMAX 
WRITE (6*1 05)  J.TH(J)*OETT(J) 

105  F ORMAT (19X*I2*6x*F6.2*9X*Fb.3) 

TH(J)«TH(J)*DTR 

104  CONTINUE 
GO  TO  41 

43  CONTINUE 

C  RE  AO  AND  WRITE  CONTROL  POINTS  FOP  NOSETIP  SMAPF 

CALL  SHAPE 

IF (IGEOM. NE. 4)  GO  TO  201 
REAP  (5*  102)  (THU)  «J>1*JMAX) 

WRITE  (6*203)  <TH(  J)  .JM.JMAx) 

203  FORMAT (2X.WREAD  IN  TH(J)  IN  DEfaPFE  FOR  INDENTFO  N0SFTIP*./.20X. 

1 10F10.6) 

DO  202  J>1*JMAX 
202  TH(J)«TH(J)»DTR 
201  CONTINUE 

WRITE  K>«304) 

304  FORMAT(*0*. ‘INDENTED  NOSETIP  ShapE*»/»?2X*«J*»11X**TMFTA*»10>« 
1*RB*.I3X*«XB«*  13X**YB»*UX. ‘DELTA*./) 

DO  61  J* 1 »  JNM 
CT»TH  (J) 

CALL  CTXY(XAA»YAA*DED»BF.CT> 

DETT (J)*DED 
CTT«CT/OTR 

WRITE (6*303)  J.CTT.BF.XAA.YAA.UFD 
303  FORMAT (20X*I3*5FI5.5) 

61  CONTINUE 


NSWC  TR  82-286 


41  CONTINUE 

BODY  POINTS  FOR  SPHERE-CONE 
WRITE  16(306) 

306  FORMAT (*0*«*COORDINAT£6  FOR  UNIT  SPHERE-CONE*,/) 

DO  58  J*1 * JNM 

X(J,l)*l.-COS(TH( J) ) 

Y(J,1)*SIN(TH(J) I 
RSH*0 • 

IF(IGEOM.EO.O)  DETT(J>*0. 

WRITE (6, 303)  J,TH(J) ,RSH,X(J,1) ,y(J.l) 

sa  continue 

DTHMIN*TH(JNM)-TH(JNM1) 

DO  56  JsJNMP.JMAX 
TH(J)*TH(J-1)*0TH 
IF(J.G£.JNM5)  TH( J) =PI? 

UETT ( J) =0ETT ( JNM) /(COS(TH( J)-TH( JNM) ) ) 

X ( Jt 1 > =DTHMlN*CQS (TMM) ♦X(J-l.l) 

y(j,i)*dthmin*sin(tmm)*y(j-i,i> 

STREAMWISE  COORDINATE  STRETCHING  ON  CONE  PORTION  FOR  J  GT.  JNM 
0TNMIN=CC  *0THMIN 

WR1 TE (6 .307)  J*  TH(J),X(J»1)»Y(J»1)»DETT(J) 

307  FORMAT (20X, I3.F 15.5. 15X.3F15. 5) 

56  CONTINUE 

SHOCK  LOCATION 
XM2=XMACH**2 

DLTOs(GAMMl*XM2*2.)/( (GAM*1.)  *xm2)*0.78 
IF ( XMACH-2.5)  21 .22.22 

21  WRITE 16.191 ) 

191  FORMAT (1 HO »*NOT  READY  FOR  M1NF  LESS  THAN  2.5*) 

(.0  TO  81 

22  IF (xmACH-10.)  23.23.24 

23  YX1*2. 376-0.  1834*XMACH*0.01036*XM2 
GO  TO  25 

24  YX1*1 .576-0.0018* (XMACH-10. ) 

25  RB=YXl**2*0.S/(l.*0LT0) 

UO  7  J=2. JMAX 

IF  ITH(J) .EO.P12)  GO  TO  71 
TCH*TAN(TH(J)) 

CGC=X ( J, 1 ) «Y ( J, 1 ) /TCB*OLTO 

Y ( J.KMAX ) *-PB/TCB*SQRT ( (PB/TCB) **2*2.*RB*CGC) 
X(J,KMAX)*X(J,1>-(Y(J,KMAX)-Y(J,1 ) >/TCH 
GO  TO  7 

71  V(J,KMAX)*S0RT(2.*P8*(X(J»1)*DLT0> ) 

X  < J.KMAX) *X (J. 1 ) 

7  CONTINUE 

X ( 1 «KMAX) =X(2»KMAX) 

Y(1,KMAX)*-Y(2.KMAX) 

WRITE  (6,701 )  .  „  .. 

701  FORMAT (/,*0START1NG  BODY  AND  BOW  SHOCK  LOCATIONS*./) 

125  FORMA TU5)U*XB*.1BX.*YB*. 18X. *XS*»lBX.*VS*.16X.*THtTA*,l4X.*J*) 

00  9  1, JMAX 

X(J,l)*X(J.l) *OMFGA 
Y(J,1)»Y(J.1)*0MFGA 
X (J.KMAX) «X ( J.KMAX >*OMEGA 

Y (J.KMAX) *Y (J.KMAX) *OMEGA  .... 

WRITE (6, 124)  X ( J. 1 ). Y ( J.1).X( J.KMAX). Y( J.KMAX), TH(J).J 
124  FORMAT (5F20.6, 15) 

9  CONTINUE 

IF (IGEOH.EO.O)  GO  TO  64 

FACTA-O.  ______  _ _  _ 

FACTA  IS  THE  FRACTION  OF  UELTA  Al  READY  DtFOBMEOt  FACTS  IS  THE  FPAC 
OF  OELTA  TO  BE  DEFORMED  IN  THIS  PUN.  FACTT*FACTA*FACTB 
READ(5«102)  FACTB 

factt»facta*factb 

WRITE (6,305)  FACTA, FACTB.FACTT 

305  F ORMAT (*0*«4X,*FRACTI ON  OF  DELTA  PREVIOUSLY  OONE«*»F5.?./»5X. 
l*FRACTION  OF  DELTA  FOR  THIS  RUN**»F 5. 2./»5X»*TOTAL  FRACTION  OF 

20ELTA  completed** .f 5 .2,/) 

DO  63  J*1,JMAX 
63  OETL(J)«OETT(J)*FACTB 


uouguu  ggoouu 
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64  CONTINUE 

3  CONTINUE 

C  FILL  ET*  COORDINATE  STRETCHING  ARRAY 

CALL  ETATB(ET.CF.KMAX) 

C  DETERMINE  X  AND  V  FOR  GRIO  POINTS  BETWEEN  BODY  AND  SHOCK 
IF ( IGEOM.EQ. 1 )  WRITE (6* 106) 

106  FORMAT </»2X»*READ  IN  XB(J> ,YB(J> ,XS(J) .YS(J)*./.20X.*J*.SX.*X8( J>* 
l,SA,*YB(J)*.5X,*XS(J)*.5X»*YS(J>*> 

DO  5  J*1 • JMAX 

C  READ  IN  XB. YB.XS. YS 

IF ( IGEOM.NE.  1)  GO  TO  4 
READ (S* 102)  XB.YB.XS.YS 
WRITE  (6, 108)  J.XB.YB.XS.YS 
108  FORMAT (1SX.I2.4F 10.3) 

GO  TO  6 

4  CONTINUE 
XS=X ( J*KMAX > 

YS=Y ( J.KMAX) 

XH=X(J.l) 

YB=Y ( J»  1 ) 

6  CONTINUE 
OXX=XS-XB 
UYY=YS-YB 
DO  5  K=1,KMAX 
ETAaET  <K ) 

X(J,K)=XB*OXX*ET* 

Y(J.K)*YB*DVY*ETA 

5  CONTINUE 
81  CONTINUE 

RETURN 

ENU 


SUBROUTINE  INI T I A 

COMMON/COM1 / JMAX  »KMAX  »  JM.KM. XMArM.GAM.OAhMl .CN.DT , SMU. JCS.PRT . 

1  iPRTtH, OMEGA, I T tTAU. I  TER. 1  NT. RThRT»R1mF, PINF.OINF.CINF.PT. ITS* 

2  1W1.IW2, I AFHD, IGEOM, TM, IVIS.IfPAN.CF »CC» JNW.RE  Y.PHD.CVlS.CVISl, 

3  TWA.ITWA.LIP.KPFS.SMJIMP.MTINE  .FTINF.SlNF.EI  INF  »  RE  YIN,  Si 'M  (40 )  » 
4DETT (40) .OETL (40) »t T 1 40) . TH C40) , T TF .F ACTE .FACTT .RE YNLO.PPTURR 

COMMON/COM2/X (40.40 ) , Y (40 ,40 ) ,XFX (40.40.2) ,XFY (40,40, P) .0(40,40) 
C0MM0N/C0M3/Q (40.40,4) ,EF (40 .4) ,c (40,40.4) ,S(4) ,A8 (4.4) .hVFC (40,4 ) 
COMMON/C0M4/A (40,4,4) ,B(40,4,4) ,0(40.4.4) .HD(4fl,4«4) • 

1UO (40.4,4) , AX (40 ) .AY (40) ,»X (40 ) ,PY (40) 

COMMON/ V I SK /CMUK AP ( 40 ) , TURMU ( 40 . 4  0 ) 

DATA  AX/40*0.0/,BX/40*1,0/.AY/40«0.0/,HV/40*) ,0/ 

DATA  TURMU/ 1600*0,0/ 

C  THIS  SUBROUTINE  INITIALIZES  THE  FLOWFIELD 
PI=4.*ATAN(1.) 

READ (5, 108)  JMAX.KMAX,ITER»IPRT,TR1.IW2»IAFBD 
108  FORMAT (815) 

JMAX*TOTAL  plints  In  j-arkay.  kmax*total  points  in  k-array 

1 TER«TOTAL  INTEGRATION  STEPS 

IPRT*1  FOR  DETAILED  PRINTOUT  IN  FlGEN.  *0  OTHERWISE 
IR1*1  READ  DATA  FROM  TAPE1.  *0  OTHERWISE 
IW2*1  WRITE  DATA  ON  TAPE 2  FOR  RESTART.  *0  OTHFRWlSE 
IAFBDsl  STORE  STARTIN3  DATA  FOR  AFTERBODY  CAl  *0  OTHFRWlSE 
READ (S, 108)  JNM, IGEOM.L IP.KRES. I TRAN. IV IS 

JNM* JUNCTURE  OF  SPHERE  ANU  CONE.  LIP*  NO.  OF  STEPS  TO  COMPLETE  THE 
DEFORMATION.  KRES=PRlNTOUT  INTERVAL  IN  K  ARRAY  FOR  RESIDUE  INFORMA 
I TRAN*NO.  OF  POINTS  OVER  WHICH  TMETA  3FCOMFS  PI/2,  MUST  HF  LT , JMAX 
IGEOMwO  UNIFORM  SPACING  ON  NOSE.  *1  READ  IN  XB.Yfl.xS.YS. 

*2  CAL.  DELTAS  AND  FINAL  XB.YP.  «3  PEAO  IN  DELTAS 
IVIS=0  INVISCIO  FLOW,  *1  LAMINAR  FLOW 
JMbJMAX-1 

km*kmax-i 

RE AO (S. 1 07)  XMACM.GAM.TM.OMEGA.CM.CF.CC.SMU.SMlllMP 
107  FORMAT (7F10.0.2F5.0) 
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c  xmach*eree  stream  mach  number 

C  GAM*  RATIO  OF  SPECIFIC  HEATS.  TM*CON» ( AFTERROD Y>  HALF-ANGLE 

C  OMEGA*  RADIUS  OF  SPHERE-CONE •  *  0  FOR  ADDING  POINTS 

C  CN* INPUT  COURANT  NO. 

C  CF*STRETCHING  PARAMETER (BtTA)FOR  GRID  POINTS  IN  K 

C  SMU*EXPLICIT  DISS.  COEF.»  SMUIMP*IMPL1CIT  DISS.  COtF. 

.  WRI TE (6. 1 0?) 

102  FORMAT (*1*»2X»*AX1SYMWETRIC  FLOW  OYER  NOSETIP*.//) 

WRITE (6.103)  XMACH,GAN.TM.0MFGA.|RlfIW2tlPRT»IAFHDtIGFnM 

103  FORMAT (*0*,2X»*MACH  NUMBER  =*»FS.2./,3X»*RATin  OF  SPECIFIC  HEAT  *• 

1 ,F5.2,/«3X,*C0NE (AFTERBODY)  HALF-ANGLE  «*»F7.3»2X»*DEGREES*./»3X. 
?*OMEGA  **»F7.3»5X,*(O*IEGA.GT.0»Om£GA  IS  THE  RADIUS  OF  SPHERE-COME* 
2  IF  IGE0M=30R4  OMEGA  VALUE  IS  RECALCULATED*, /,?IX,*IN  SUP.  SHAPE* 
30MEGA*O.MORE  RAYS  TO  BE  ADDED) * ,//.3X»*IRl  **.I2.5A.*(  1  FOR  READ 
4TAPF1S  0  OTHERWISE)*. /.3X»*|W2  **.I2.SX»*(1  FOR  WPTTF  ON  TAPE? 
5*  0  OTHERWISE)*.  /.3X.*IPRT  **,I2.SX.*(  1  FOP  DETAILED  WRITE  OUT 
6  FROM  EIGEN*  0  OTHERWISE ) *,/»3X»*lAFbO  **»I?.SX.*(  1  FOR 

7ST0RAGE  OF  STARTING  DATA  FOR  AFTFRBOOY  CAL.*  0  OTHERWISE)** 
8/.3X.*IGE0M  **» I2.5X  »* (  0  FOR  UNIFORM  SPACING  ON  NOSE  *  )  FOR  RFAD 
9  IN  XB.YB.XS.YS  S  2  FOR  READ  IN  THCJ)  AND  DETT(J)  **»/.17X. 

1*3  FOR  CAL.  DELTAS  AND  FINAL  xH.Yd  WITH  UNIFORM  TH(J)*  4  FOP  REA 
20  IN  THU)  AND  CAL.  FINAL  XR.Yb)*) 

WRITE (6.207)  LIP. 1 VIS. CF.CC.l TRAM. KRES.Smii.SmUIMP.cn 

207  F  ORMAT (*0*»  2X.*LIP  **.I4.5X.*(  0  FOR  WITHOUT  SHAPE  CHANGE  * 

3n  for  shape  change  completed  in  n  stfps)*./.3x.*ivis  **.i?.*»x.*( 

40  FOR  INVISCID  FLOW  S  1  FOR  LAMTNAR  FLOW) • .//*3X »*CF (BETA) *• .FI ?.S 
S.SX.*(  FOR  UNIFORM  SPACING  SET  TO  10000)*. /.3X.*CC  **.F5,?.SX.*( 
6STRFTCHING  FOR  POINTS  BT.  JNM* I TDAN  AND  JMAX » *./♦ 3X.*I TPAN  **.IP» 
75X»* (MUST  HE  LT.JMAx-JNM  FOP  THFTA  TO  GO  TO  PI/?) *./.3X.*KRES  **» 
BI2.5X.*( INTERVAL  IN  K  FOP  RESIDMF  INFORMATION) *»/»3X .*EXPLICIT  DIS 
9SI.  COEF.  **«F5.3./. 3X.*IMPL ICi T  DISSI.  COEF.  **.FS.3*/.3X, 

1 *COURANT  NO.  **.F6.2.//> 

WRITE  (6.208)  JMAX. KMAXUNM.  ITER 

208  FORMAT (*0*»2X»* JMAX**, I5./.3X  «*WMAXS*. IS./.3X ,*JNN**. IS.SX. 

1  ‘(JUNCTURE  OE  SPHERE  ANO  CONE) *,/.3X.*I TER  =*, I4.SX.* (TIME  STEPS 
2F  OR  THIS  RUN)*) 

GAMM1=GAM-1 
GAH1I*1.0/GAMM1 
I TF  =0 
TAU=0. 

IT*0 
ITS=1 
F  ACTT*0. 

FACTA*0. 

FACTH*0. 

JCS*1 
PINF*1 • 

R I NF* 1 . 

CINE  *SORT (PlNF*GAM/Rl Nf > 

UlNF*XMACH*CINF 

IF ( IVIS.EO.O)  GO  TO  9 

C  RE AO  CONSTANTS  NEEDED  FOR  VISCOUS  FLOW  CAL. 

READ (5* 104)  REY.PRD.PRT.CVIS.TWA.ITWA.ITUR.ITE 

104  FOHMAT(5E10. 0,315) 

C  REV*EREE  STREAM  REYNOLDS  NO.,  PPD*  FREE  STREAM  PRANOTL  NO. 

C  PRT*  TURBULENT  PRANOTL  NO.,  CVlS*CONSTANT  IN  SOUTHtRLAND"S  LAW 
C  FOR  VISCOSITY.  TWA*  WALL  TEMPERATURE 
C  I TwA*l  ISOTHERMAL  WALL,  *0  ADIABATIC  WALL 
C  I TUR*1  TURBULENT  FLOW,  *0  LAMINAR  FLOW 

C  ITF*1  PRINTOUT  STANTON  NO.  ONLY,  *2  PRINTOUT  T -FIELD  ALSO**0  NO  H 
REYIN*REY 
REYNLD*REY/OINF 
PRTURb*PRD/PRT 
CVIS1*CVIS»1. 

WRITE (6, 105)  REY.PRD.PRT.CVlStlTWA.ITUR.ITF 

105  FORMAT (*0*.?X.*RE  **,E15.6,/,3X,*PR  «*.Fb.3./.3X.*PR(TllRH.)  «*.E 
16.3, /,3x,*C VIS  *  110/TINF (KELVIN)  **.FP. J.5X ,* (  CONSTANT  USED  IN  S 
2UTMERLANOMS  LAW  OF  VISCOSITY  >*,/,3X,*ITwA  *•♦!?, 5X,*<  0  FOR  AD TAB 
3ATIC  WALL*  1  FOR  ISOTHFRMAL  WAL|.)*»/*3X,*ITUP  **.I2.5X,*(  0  E OR  LA 
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130 


135 


140 


145 


150 


155 


160 


165 


4m1NAR  *  1  FOR  TURBULENT )*./.3X«*TTF  «*.I*:.SX.*<  1  FOR  PRINT  OUT  5T 

5  NO.  ONLY  S  2  PRINT  OJT  T-FIELO  ALSO)**//) 

IF ( ITWA.EQ. 1 )  WRI Tt (6* 1 06)  TWA 

106  FORMAT (*0*. 2* «*TI»  =*FB.3.//> 

9  CONTINUE 

C  SET  UP  CONSTANTS  AT  FREE  STREAM 
WRITE  16*109) 

109  FORMAT <*0*.«FREE  STRt AM  CONDITIONS*) 

UlNFsQlNF 

VlNFsO. 

HT I NF WGAM/GAMM 1 »P INF/ R I NF ♦ 0 . 5*0 1 mF  **2 
fc  TINF*HTINF-PINF/PINF 
SlNF  sPlNF/P INF  **6AM 
1 1 1 NF  * 1 . /GAMMl *P I NF /K 1 NF 

WRITE <6(100)  PINF (RlNF(QINF(ClNF .UINF  * V INF ( HTInF (ET INF *S t mF.E I  IMF 
100  FORMAT (*0*.2X.*PInF (PRESSURE)  «*. Fh. 4. /.j*.*PINF (DENSITY)  =*.F«.4, 
2/*3x**OINF (TOTAL  VEL.)  »*.F 8.4,/, 3X.*A1NF (SOUND  SPtED)  **.F8.4./. 
33 A, *U INF (U  COMP.)  **.FH.4./«3X,«vlNF (V  COMP.)  =*.Fb.4./. 

4  3X(*mtinf<t.  enthalpy)  =*.fb,4./.3x.*etinf (T.  spec,  fn 

4ER0V)  **»F8»4,/»3X .*SINF (tNTPOPY)  **.F8»4./»3X.*£I INF ( INTFRNAL  FNE 
6RGV)  **«Ffl.4»//) 

CALL  ETATB(ET.CF.KMAX> 

WRITE <6«1 12) 

WRITE <6. Ill) <ET(R)(K>1.KMAX) 

112  FORMAT (*0*.2X, ‘NORMALIZED  DISTANCE  FROM  BODY  TO  SHOCK*) 

111  FOHMAT(?0X(10F10.6> 

Xl=(2.0*GAM*XMACH«*2-SAMMl)/(GAM*1.0) 

X2- (GAM* 1.0) *XMACH**2/ (GAMMl *XMACH**2*2. 0 ) 

Pl=Xl*PINF 

Hl=X2*PINF 

ENT*Pl/Rl**GAM 

PT= ( 1 ,0/X  t ) *• ( 1 .0/GAMH1 ) * (0.5*  <G»M* 1 . ) *XMACH*«?) ** (GAM/GAMM1 ) *PTNF 
XX=1 .0*0.5*GAMM1*XMACH**2 
PTORT»XX*PINF/HINF 
WRITE (6.117)  PT 

117  FORMAT (/(2X#*STAGNATI0N  PRFSSURF  PT«**E10.4) 

C  CHECK  FOR  FRESH  START  OR  CONTINUATION 
IF (IRl.EQ.l)  GO  TO  22 
CALL  GRID 
CALL  JACOB 

C... INITIALIZE  0  VECTOR  TO  FREE  STREAM  VALUES 
DO  1  KM.KMAX 
DO  1  J= 1 • JMAX 
01*1 .0/0 (J.K) 

0<J(K(1)=RINF*OI 

0(J.K(2)sRINF*UINF*DI 

0<J,K.3)=RINF*VINF*DI 

0(J.K(4)«(PINF*GAM1I*RINF*OINF*»?*0.5)*DI 
C...SET  S  ARRAY  TO  0  EVERYWHERE 
DO  1  N»1 .4 
1  S ( J.K (N) =0 . 0 

C... INITIALIZE  FLOW  FIELD  FOR  BLUNT  BODY  PROBLEM 
GAMP1*GAM*1.0 
DO  ?  J*2. JMAX 

IF  (ABS(XEX(J. 1.2)) -0.000001)  6(6*7 

6  TMET*0i5*PI 
00  TO  8 

7  THET=ATAN  (XEY( Jt 1  (?) /XEX ( J. 1 *2) ) 

6  CONTINUE 

K*KmAX 

SANC,=0.5*PI-ATAN(-XFV(  J«K*2)/XEX  (J.K. 2) ) 

XX»XMACH**2*5IN(SANG) •*? 

PS=  <2.0*GAM*XX-GAMM1 ) /GAMP 1*P INF 
RS=GAMP1*XX/ (GAMMl*XX*2.0) *R1NF 
US»(1.0-2.0*(XX-l.0)/3AMPl/XMACH**2)*3INF 
VS«?. 0*(XX-1.0)*C0S (SANG) / (GAMP! *XMACH**2*SIN (SANG) )*OINF 
PH*PINF*( (PT/PINF -1.0)* (1.0-1.0?*SIN( ThF I )•*?♦(). 12*SIN(THFT)**4)* 
W  1.) 


170 
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RB*  (PB/ENT ) ** ( 1  •  O/GAM) 

'  OB«S<)RT(2.0*GAM/GAMMl*ABS(PTOKT-PH/kM) >  - 

YY«PI*0.5-THET 
UB»AHS (OB*COS  < Y Y) ) 

VH«0B*S1N(YY) 

DO  ?  K*I «KMAX 
YV*FT<Ki 

PHtSS»PB*YY*  (PS-PH) 
kHO«RB*YY*(RS-RB> 

OVE LN*SQRT  ( 2 . 0*GAM/GAMM 1 *Atf S  < PTOPT -PRESS/RHO) ) 

UVEL*UB*YY*(US-UH) 

VVEL=VB*YY*(VS-VB> 

QVELOOQRT  <UVEL**2*VVEL*«2> 

PATOVELN/OVELO 

uvel=uvel*rat 

VVtL*VVEL«kAT 
1)1  =  1  *0/0  ( J*K ) 

Q(J,K.l)*RHO*DI 

0(J,K«2)=RH0*UVEL*UI 

G(J.K»3)sRHO*VVEL*0I 

2  O(J.K»6)s(PRESS*GAM|I*RH0*(UVEL**2*VVEl  **2>*0.5>*DI 
c... reflect  metrics  ano  DEPENDENT  VARIABLES  about  plane  OF  SYMMETRY 
00  6  MlfKMAX 
DO  ,K) *0 (2»K) 

XtXU.K,l)*-XFX(2#K,l> 

XEY(l*K»n»XEYC2.K.l) 

XEX(ltK*2)>XEX(2*Kt2) 

XEY ( 1 *K*2I s-XEY (2.K.2) 

DO  S  N*1 *4 
S  (><1,K,N)*Q(2.K,N> 

A  0(l(Ki3) *~Q (2*K  *3) 

GO  TO  26 

22  CONTINUE 
REWIND  1 

REAOO)  (  (XlJtKI  »  J*1  .  JMAX)  »K«1  «KmAX)  • 

1  ( IV(JiK) »J*1.JMAX) »K«1»KMAX) • 

1  ((D(J«K)»J=1» JMAX) »K*I »KMAX) * 

I  <( <XEX( J.K.N) .J-l.JMAX) ,K*1.KWAX> ,N=1.?>  t 

1  (( I XEY (J.K*N>  *J*1* JMAX) ,K*1*KHAX>  ,N«1,2> , 

1  ( ( (0 ( J*K*N) . J«1 1 JMAX >  *K*1 »KMAX> ,Nsl ,4) , 

1  JMAX.KMAX.XMACh.GAM* IT »TAU*F  ACTA*  (DFTT ( J) tJ«l*JHAX) 

XMACH«0 INF/C INF 
ITS*IT*I 
ITER«ITER*IT 
WRlTE<6,110) 

110  F ORMAT I  *0*  »  *5TART I N3  SOLUTION  WAS  READ'  FROM  TAPE*) 

C  CHECK  FOR  OPTION  OF  ADOlNG  POINTS 
IF (0ME6A. EQ.O)  GO  TO  23 
C  CHECK  IF  FURTHER  DEFORMATION  IS  NEEDED 
IF(LIP.EQ.O)  GO  TO  21 

C  FACTA  IS  THE  FRACTION  OF  DELTA  A|  READY  DEFORMED.  facTB  IS  THE  FRAC 

C  OF  OELTA  TO  BE  DEFORM  ED  IN  THIS  PUN*  rACTT«e'ACTA*rACTB 

REAO(S*107)  FACTB 
FACTT«FACTA*FACTB 
WRITE<6*305>  FACTA, FACTB. FACTT 

305  FORMAT(*0**6X,*FRACTIJN  OF  DELTA  PREVIOJSIY  OONF»*,!-S.^,/,«jX, 
I*FRACTION  OF  OELTA  FOR  THIS  RJN-#.FS.?./,SX.*T()TA_  FRACTION  OF 

2delta  completed***^  5.2*/) 

CALL  JACOB 
DO  306  J* I .JMAX 

306  OETUJ)*OETT(J)*FACTB 
21  GO  TO  26 

C  OMEGA  *  0  ADDING  GRID  POINTS 

23  CONTINUE 
CALL  1NTFRP 

26  CONTINUE 

SUM (2) "SORT (X(2»1)**2*Y(2«I) **2) 

DO  II  JO* JMAX 

II  SUM(J»«SUM(J-11*SQRT((X(J,1)-X(  i-1 .1) ) **2* (Y ( J,1 ) -Y< J-l .1) )••?) 
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30 


35 


WRITE (6*113) 

113  FORMAT (*0*«»ARC  LENGTH* ) 

WHITE  <6*1 14)  (SUM(J> *J«2.JMAX) 

114  FORMAT (20X»10F10*5) 

WRITE  16*401 1 

401  FORMAT (/»*0STAHTING  FLOWFIELD  INFORMATION**/) 
RETURN 
ENO 


SUBROUTINE  INTEGR 

COMmON/COM1/JMAX*KMAX» JM*KM«XMACH*GAM»OAMMl »CN»DT  *SMU» JC5*PRT  * 

1  IPRT»H»OMfcGA»IT»TAU*ITER*ENT*PTCRT «P1NF*RINF  *QINF»CINF*PT* ITS* 

2  1R1 ♦ IW2* IAFBD* IGEOM*  TM*IVIS*IlRAN*CF *CC* JNM*RFY*PRD*CvIS*CVIS1 • 

3  TWA«ITWA*LIP*KRFS*SMJIMp*HTINF *FTIN) . SINF .El INF, HEVIN.SUM (40 ) * 
40ETTI40) tOETL (40) *ET(40>  »TH(40) *lTF*FACTb*FACTT *REYNLO*PPTURb 

C0MM0N/C0M2/X (40*40) *Y (40*40) *XFV (40*40*2) *XFY (40*40*2) *0(40*40) 
COMMON/COM3/O(40*40*4) *EF (40*4) (40 *40 *4) *6(4) * AS (4*4) *HVEC (40*4) 
C... COMPUTE  FORCING  FUNCTION  AND  STORF  TEMPORARILY  IN  S  ARRAY 
CALL  RHS 

C... COMPUTE  RESIDUE  EVERY  25  STEPS  TO  CHECK  FOR  CONVERGENCE 
IF(MOD(IT*2S).EQ.O)CALL  RESIDU 
C...AOO  FOURTH  ORDER  DISSIPATION  TO  SMOOTH  SOLUTION 
CALL  DISSIP 

C... SOLVE  FOR  O-BAR-BAR 
DO  1  K*2*KM 

C...FILL  ELEMENTS  OF  I*H*DX  A  FOR  BLOCK  TRIDIAGONAL  INVERSION  AT  EACH 
C...K  TH  LEVEL 

CALL  LBLTRA(K) 

C... INVERT  BLOCK  TRIDIAGONAL  MATRIX  AT  K  TH  LEVEL  AND  STORE  SOLUTION  IN 
C...S  ARRAY 

CALL  BTRI (2* JM) 

DO  1  L«l«4 
DO  1  J*2«JM 

1  S( J*K*L) *EF ( J*L) 

C... SOLVE  FOR  O-BAR 

DO  2  J*2* JM 

C. . .FILL  ELEMENTS  OF  I*H*DV  B  FOR  BLOCK  TRIDIAGONAL  INVERSION  AT  EACH 
C...J  TH  LEVEL 

CALL  LBLTRB(J) 

C... INVERT  BLOCK  TRIDIAGONAL  MATRIX  AT  J  TH  LEVEL  AND  STORE  SOlUTION  IN 
C...Q  ARRAY 

CALL  BTRI (2*KM) 

DO  2  L«1.4 
00  2  K>2»KM 

2  0( J*K*L) *EF (K*L) *Q( J»<»L> 

TAU»TAU*DT 

RETURN 

ENO 


10 


15 


SUBROUTINE  INTERP 

COMMON/COM 1 / JM AX  *  KMAX *  JM  »KM  «  X MACH  »  GAM*  G AMM 1  *  CN  *  DT  * SMU*  JCS • PRT • 

1  IPRT*H*0MEGA.IT»TAU* I  TER* ENT  »PTORT  *PINF»PINF *01 NF  *CINF  »PT» ITS* 

2  1R1 ♦ IW2* I AFUD* IGEOM*  TM* I VIS* I TRANtCF  »CC» JNM»RFV»PRD*CVlS.CVlSl • 

3  TWA»ITWA*LIP*KRES*SMJIMP.HTINF,FTINF.SINF*EIINF.REY1N,SUM(40)* 
4UETT (40)VDETL (40) *ET (40) *TH(40) *TTF.FACTB.FACTT*REYNL0«PRTUPB 

C0MM0N/C0M2/X (40*40) »Y (40*40) »XFX (40*40*2) *XEY (40*40*2) *0 (40*40) 
C0MM0N/C0M3/Q (40*40*4) »EF (40*4) (40*40*4) *G (4) * AS (4*4) *HVEC(40*4) 
COMMON/BDTH/XOl *X02*X03*X04*Y01 *y02*Y03« Y04*SL1 »SL2*SL3.P1*R2» 
1R3*R4*CT1 »CT2»CT3*CT4*CT5.CT6»X00»RS0DV 
DIMENSION  P(40) * YA (20) *XA(20) *THAD(20) 

DIMENSION  XZ (40) * YZ (*0) *QZ(40*4) ,DZ (40) *ETZ(40) 

THIS  SUBROUTINE  INTERPOLATES  FLOW  VARIABLES  FOR  NEW  GRID  POINTS 
PI>ATAN(1.)*4. 

DTRaPl/180* 

P2»0.5*PI 
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RE ADIS. 100)  JAA. JIX.KIM 

100  FORMAT  (3IS) 

IF(JAA.EQ.O)  60  TO  51 

20  WRITE (6«138) 

138  FORMATI*0*.*INFORMATION  FOR  NOSFTIP  SHAPE*) 

WR1TE16.200)  JAA.JIX 

200  FORMAT t*0*.*ADDING  GRIDS  IN  J-ARRAY.  NO.  OF  RAYS  ***I2.lOXt*JIX  •• 

1.12) 

25  C  RE AO  INPUT  VALUES  OF  JAA  THETAS  TO  BE  ADDED 

READ (5. 101 )  (THADIJ) »J*1*JAA) 

101  FORMAT (8F10.0) 

C  READ  AND  WRITE  CONTROL  POINTS  FOR  NOSET1P  SHAPE 

CALL  SHAPE 

30  WRITE  16.205) 

205  FORMAT <//.*0NEW  POINTS  ON  BOOYt*) 

00  9  N*1 .JAA 
CT*THADIN> 

CALL  CTXY(XAA.YAA.OED.BF.CT) 

35  XAIN)*XAA 

YA(N)*YAA 
CT*CT/0TR 

WR1TE16.201)  N*CT.XA(N).YA(N) 

201  F0RMATI5X»*RAY  ■*.I2.5X.*AT  THB  ***F7.?*SX«*XA  ■**F8.4«SX.*YA  »*.F 

40  18.4) 

9  CONTINUE 
JS*2 

00  495  J-2.JMAX 

495  TH(J>*-ATANI(Y(J»KMAX)-YIJ.l))/<X<J»KMAX)-X<J.l) J ) 

45  THIDW-THI2) 

DO  11  1*1 .JAA 
JMAX*JMAX«1 
JNMaJNM.l 
00  12  JaJS.JMAX 

50  IF(TH  (J).LT.THAO(D)  GO  TO  12 

JbFaJ-1 
JAf*>l 

RATA*ITHIJ>-THAD<I>)/ITHIJ)-TMlJRF)l 
RATE* ( TH U) -THAO II ) ) / 1 TH (  JAF ) -TM I J)  ) 

55  HO  3  JAsJAF.JMAX 

J I * JMAX- JA* JAF 
JL*JI-1 
DO  3  K*1 .KMAX 
X IJI *K)*X I JL.K) 

60  Y IJI *K) *Y (JL.K) 

0«J!.K>*01JL.K) 

TH(JI)«TH(JL> 

00  4  N*l«4 

4  0(J1.K.N)*Q<JL.K.N) 

65  3  CONTINUE 

TMI J)*THAOI I ) 

THIi»*-THI2) 

OA*SORT(  IX(J-1»KMAX)-XIJ-1*1) )**?MY<J-l*KMAX)-YlJ-l*l>>**2> 
OH*SORTI (XI JM.KMAX) -XI J»l.l ))•*?♦! Y I J*1 .KMAXI-VI JM .1 > >**2> 

70  US*©84RAT  A* IOA-DB ) 

5X*0S*C0S  (THU) ) 

SY*OS*SINITHI J) ) 

00  5  K*1 .KMAX 
XlJ.K)*XAtI)-SX*ETIK) 

75  YIJ,K)*YAII)»SY*ETIK) 

JFF-JAF.1 

PJ*0(J.K.I )*DIJ*K) 

RJbaQ I J0F  »K. 1 ) *0 ( JBF •<) 

RJA*OIJFF.K.1)*OIJFF.IO 
80  UJ*QI J.K.2) /Q( J.K.l ) 

l)J0*OI  JBF  .K.  2  I/O  I  JBF  .X.  1 ) 

UJA*0IJFF.K.2)/0I JFF.X.l) 

VJ*0( J.K.3)/0I J.K.l ) 

VJ0*O1 JPF.K*3)/0IJBF*<.1) 

85  VJA*01JFF.K.3)/0(JFF.<*1) 

E J*0( J.K.4)*0I J.K) 
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p: 


ft 


90 


95 


100 


105 


110 


115 


120 


125 


130 


135 


140 


145 


150 


155 


EJUsQi JHF .K.4)  *1)1  JbF.K) 

EJAsQ<JFF.K.4)*D< JFF.O 
IF(I.GE.JIX)  GO  TO  15 
RJ=RJ-RATA*1RJ-RJB> 

U J«U J-RAT A* ( U J-U JB 1 
VJ=VJ-RATA*<VJ-VJb) 

EJ=E J-RATA* IE  J-E  JB) 

GO  TO  16 

15  rj=rj-rate*irja-rj> 

UJ=UJ-RATE*  1UJA-UJ) 

VJ=VJ-RAT£*1VJA-VJ> 

EJ=EJ-RATE*1£JA-EJ) 

lb  0 1 J.K) =D I J.K) -RATA* (D 1 J.K) -D ( JHF .K ) ) 

0I=1./0<J.K) 

0( J*K.l)sRJ*DI 
O ( J.K.2) =UJ*Q1 J.K* 1 ) 

U 1 J.K .3) = V J*U 1  J.K  « 1 ) 

01J.K.4)*EJ*DI 
5  CONTINUE 
JS=J 
GO  TO  11 
12  CONTINUE 
11  CONTINUE 
JMsjMAX-1 
00  17  K=1*KMAX 
00  18  N=1.4 
18  0(1,K.N)=0(2.K,N) 

0 ( 1 .K .3) s-0 ( 2  .K .  3 ) 
l>  1 1  «K)  *U  I2.K) 

X<1,K)=X12.K) 

17  YlliK(s-Y(2.KI 
00  301  J= 1 .  JMA A 
00  301  K=1,KMAX 
00  301  N=1.4 

301  QlJ,K.N)*QlJ.K.N>*OtJ.K> 

CAUL  JACOB 

00  302  Jsl.JMAX 
DO  302  K*1.KMAX 
00  302  N=1.4 

302  01J,K»N)=U1 J»K»N)/Dl J.K) 

K*1 

00  8  J*1 *  JMAX 
Z*1./01J,K,1> 

Rsftl J.K. 1 ) *01 J.K) 

U*0( J*K  » 2) *Z 
V*U1J.K,3)*Z 
E*01J,K,4>*01J.K) 

8  PI J)«1E-0.5*R*1U**2.V»*2) )*GAMMi 
WRITE  16* 1 36) 

136  FORMAT (*0*»*SUHFACt  PRESSURE  0UTR1  3UH0N  Ar.TER  A 3U I NG  »MNTS*) 
WRITE16. 122) (P(J) .J*1.JMAX) 

122  FORMAT I20X.10F10.5) 

RETURN 
51  CONTINUE 

C  ADD  OR/ANO  REARRANGE  5RID  POINTS  IN  K-ARRAY  WITH  NEW  VALUE  CF1 
REAO(5.101)  CF1 
WR1TE16.202)  KIM.  CF1 

202  FORMAT  1*0*. *NEW  GRIDS  IN  K-ARRAy.. NO.  OF  POINTS  **»I3.5X,*NEw  STRE 
ITCHING  COEF.  ■*,F10,4F.5X,*NORMM  IZEO  OlSTANCF  FR3M  BODY  TO  SHOCK* 
2) 

CALL  ETAT81ET.CF.KHAX) 

00  54  t«l.KMAX 
54  ETZ 1 1 ) »ET 1 1 ) 

CALL  ETAT81ET.CF1.KIM) 

WR1TE16.203)  1ET1K) .K*1.KIM> 

WRITE  16.204) 

204  FORMAT I5X.*THE  OLD  VA.UES  ARE*) 

WR1TE16.203)  1ETZ1K) .<*1,KMAX> 

203  FORMAT II OX .1 OF 10. 4) 

00  52  J>1.JMAX 
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DO  S3  Kb1,kMAX 
XZ(K)*X(J*K) 

YZ4K)*Y<J,K> 

160  DZ(K)*0(J*K) 

DO  S3  N*l,4 
53  OZ (K.N) *04 J.K.N) 

00  55  K*2* KIM 
00  56  H»2,KMAX 

165  56  IF(ETZ(M).G£.ET(K))  03  TO  57 

57  RAT£«4£rZ4M)-Er(K>)/(£T7(M)-£r2rM-l)) 
X(J*K)«XZ(M)-RATE*(XZ(M)-XZ(M-ln 
V(J.Kl*YZ4M>-RAT£*4Y/(M)-Y7(M-ln 
l>(J,K)*OZ4M)-RATE*43Z(M)-DZ<M-lt ) 

1Z0  rmboz<h,i)*ozim) 

RMUQZ(M-1,1)*0Z4M>1) 

UM*0Z4M.?)/QZ(M«1) 

UM1*GZ  4M- 1  *  2) /OZ 4*4-1 • 1 ) 

VM*QZ  4M«3) /QZ 4M»  1 1 

175  VM1*QZ4M-1.3>/QZ4*4-1.1) 

EM*QZ4M,4)*DZ4M) 
fcHl*QZ4M-l»4)*0Z  4*4-1) 

H*  (RM-RATE* 4RH-HM1 ) ) 

U»  40H-RATE*  4UM-UM1) > 

ISO  V*  4 VM-RATE*  4 VM-VH1 ) ) 

f  *  4EH-RATE*  4EM-EM1 ) ) 

04J.K»1)*R/D4J»K) 

O4J,K.2)*0*Q(J*K«l) 

4)<J.Kf3)*V*Q4J*K»l) 

105  U(J«K«4)*E/0( J»K) 

PI* (E-0.5*R*  4U**2*V**2) ) •GAMHl 
55  CONTINUE 
52  CONTINUE 
KMAXsKIN 

190  KM*KMAX-1 

call  JACOB 

121  FORMAT 4 12F10.5) 

RETURN 

ENO 

1  SUBROUTINE  JACOB 

COMMON/COM 1/JMAX»5NAX» JM.KM,X*4ArH.bA*4,&AMM|  ,rN.0T*S*4U.  JC5.P9T* 

1  IPRT.H,OMEGAtIT»TAJ»ITERttNT*RTOST.RlMF,PI^F.OlNr'»ClNF.PTtnS. 

2  IRl*lM2*IAFb0.IGE0M,TM«IVIS.I  l<  AN.CF.CC.JN*t.PFY.*iO*CVlc.CVISl. 

5  3  T*A.ITbA.LIP.KRES.SMJIMP*HTINF.FTlNF.SlNF.mNF.TEVlN,SUM(40>» 

4UETTI40) .OETL440).ET44O).TH440>.TTF.F*CTb.FACTT.REYNLn.BBTU»B 
COMMON/COM2/XI40.40) *7(40*40) *Xrx(40» 40*2) •  *FY 440 *40*2) *0(40*40) 
COMMON/CON3/Q(40*40*4)*EF(40«4>.C (40*40*4)* 3(4)  *A-M4*  4)  *HVEC(40.4) 

oata  IFLAG/O/ 

10  1F4LIP.EQ.0)  GO  TO  13 

IF(IFLAG.EQ*1)  go  to  u 
00  26  J*1 * JMAX 

26  UETl(J)*OETL<J)/4FLOAT(LIP)) 

*RITE46«101) 

15  101  FORMAT 4*0800Y  SHAPE  CHANGE  BE  I No  INSTITUTED*) 

11  CONTINUE 

IF4ITrITS.GE.LIP. OH. IT.LT. ITS)  r.n  TD  13 
DO  14  J*1.JMAX 

IF (A8S4X4  J*l)-X(  J*K*4AX) ) .LT.1.0r-6)  G3  TO  15 
20  SLP* (Y  4 J* 1 )-Y ( J.KMAX) ) /(X  4 J* 1 ) - v  4 J»<MAX) ) 

DD*5QRT  4  SLP**2* 1.0) 

GO  TO  15 

16  CC>0. 

DO— 1.0 

25  GO  TO  17 

15  CC*1./0D 
0D»SLP/00 

17  CO)  'INUE 

X  4  ,,*X4J»'  -C*OE TL4 J) 

30  Y4J*i.*  I J,  »00*DETL 4 J) 

00  K  '  -2 «KM 
ETA*ET4X) 

X(J*Kk>(X(J«KMAX)-X(J«l))*£TA»Af.l*l> 
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35 


40 


45 


SO 


55 


60 


65 


70 


75 


80 


85 


14  V  ( J.K)« ( Y  ( J.KMAX)  -Y(Jtl))*ETA*Y(.t*l) 

13  CONTINUE 
JMMsJM-1 
KHMsKM-1 

C... COMPUTE  X-XI  and  Y-XI*  DXI  AND  OtTA  3  1 
DO  1  K*1»KMAX 
00  2  J*2t JM 

XEY(J*Kt2)*(X(J*l.K)-X(J-l.K) )*0.5 

2  XEX(J.K.2)=(Y<J*l.K)-Y(J-ltK) >*0.5 
XEY(1*K,2)*(-3.0*X(1  tK)*4.0*X(2.K)-X(3.K>  )«0.5 

XEY ( JMAX »K»2) = (3.0*X ( JMAX »K  > -4. 0«X ( JMtK ) *X ( JMM.K ) ) *0.5 
XEX(1,K.2)*(-3.0*Y(1 t<)*4.0*Y<2,K)-Y(3.K)  )*0.5 
1  XEX ( JMAX tK .2)  =  (3. 0*Y ( JMAX  tK ) -4.0*Y(JM»K)*Y(  JMM*K ) ) *0.5 
C... COMPUTE  X-ETA  AND  Y-ETA 
00  3  J= 1 • JMAX 
00  4  K=2*KH 

XEY (J»K*l)=(X(JtK*l)”X(J»K-l ) )  *0.5 

4  XEX(J*K*1)=(Y(J*K*1)”Y(J«K-1))*0.5 
XEY(J*1,1)*(-3.0*X( Jtl)*4.0*X(J*?)-X( J.3) )*0.5 

XEY (JtKMAX, 1) * (3.0*X< J.KMAX > -4. 0*X ( Jt <M> *X (J.KMM) > *0.5 
XEX ( J« 1 • 1 ) * (-3.0*Y ( J* 1 ) *4.0*Y ( J«?) -Y ( J«3) ) *0.5 

3  XEX ( J,KMAX,1)=(3.0*Y( J.KMAX ) -4. 0*Y ( J. <M) *Y ( J.KMM) )*0.5 
C.. .COMPUTE  XI-X.  XI-Y.  ETA-Xf  AND  ETA-Y 

DO  5  K*1.KMAX 
00  5  J=1 • JMAX 

DI*1.0/(XEXU.K.l)*XEY(J.K.2)-XFY(J.K.l)*XF X ( J.K.2)  1 
011=01 

IF (IFLAG.EO.O)  60  TO  7 

C.. .ADJUST  CONSERVATIVE  VARIABLES  BASED  ON  NEW  MESH 
00  6  N=I*4 

6  0(J,K.N)=Q(J.K.N)*0(J,K)/UII 

7  CONTINUE 

C...THE  GEOMETRIC  JACOBIAN  IS  DEFINED  HtKt  AND  STORED  IN  THE  0  ARRAY 
0 ( J.K) =DI I 

XEX(J.K«1)=XEX(J.K«1)*DI 

XEY(J.K,1>*-XEY(J»K*1)*DI 

XEK<J»K»2>*-XEX<J»*..2>*ni 

5  XEY(J.K.2)*XEY(J.K.2)*DI 

C... REFLECT  METRICS  ANO  DEPENDENT  VARIABLES  ABOUT  PLANE  OF  SYMMETRY 
IE ( IELAG.E0.01  GO  TO  B 
1)0  9  K=1.KMAX 
l)(l.K)*0(2.K) 

XEX(1»K»1)*-XEX(2*X*1) 

XEYI1»K»1)*XEYI2«R«1) 

XEX ( 1 »K»2)*XEX (2»K»2) 

Xf Y(1.K.2)=-XEY(2.K.2) 

DO  10  N*1  .4 
10  O(l.K«N)=0(2«K.Nl 
9  0(1»K»3)*-Q(2*K*3) 

8  CONTINUE 
IFLAG*1 
RE  TURN 
END 


1  SUBROUTINE  LBLTRA  CM 

COMMON/COMl/JMAX.KMAXtJM,KM»XMACH«GAM.GAMMl.CN.DT.SMU.JCS.PRT» 

1  IPRT tH .OMEGA* I T.TAU. ITER.ENT  .PT0KT.P1NF.RINF.QINF.CINF.PT .ITS. 

2  IK1 » IB2. I AFBD  * IGEOM» TM.IVIS.ITRAN.CF *CC* JNM.PF V.PKD.CVlS.CVISl * 

5  3  T*A.ITWA»LIP.KRES«SMJIMP,HTINF  .FTINF. SINE. EIINF. KEVIN, SllM(40>  . 

40ETT (40) *DE TL(40> »t T ( 40 >  »Th (40 > , T TF ,F ACTb.E ACTTaREYNLD.PRTURB 
COMMON/COM2/X (40.40) »Y (40»40) ,XFX (40.40.2) « XEY (40,40,2) ,0(40.40) 
COMMON/COM3/Q (40.40.4) *FF(40«4) ,9(40,40.4) »G(4) * AB (4.4) .HVECI40.4) 
COMMON/COM4/A  (40*4*4)  .8(40,4,4)  ,C  <40 , 4 .4 )' .HD  (40 ,4.4)  , 

10  lUD (40.4*4) »AX (40) »AY (40) ,BX<40) ,RY(40) 

00  1  J* It JMAX 

C. . .LOAD  BLOCK  A  MATRIX  EVALUATED  AT  N  TH  LEVEL  FOR  ALL  J  INTO  HD  ARRAY 
CALL  ABMATX ( J  K.l) 

DO  1  M* 1*4 
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15  00  1  L* 1*4 

1  HD ( J,L  *M) =AB (L  *M) 

C.  *  r  Fit  L  OFF-DIAGONAL  AND  DIAGONAL  ELFMENTS  BASED  ON  A  2-NO  ORDER 

C.. .CENTRAL  difference 

DO  2  J*2,JM 

20  SHlsSMUIHP»0(J-l.K)/0(J.K> 

SP1=SMUIMP«D(J*1.K)/0<J,K) 

DO  2  M*1 «4 
00  3  Lx 1 *4 

A(J*L*M) *-H0 ( J-l «L  «M) *H 
25  H ( J.L.M) *0.0 

3  C(J,L*M)aHD(J«-l,L*M)*H 

C...SET  8  ON  THE  DIAGONAL  REPRESENTING  THE  IDENTITY  MATRIX  TO  ONE 
A(J,M.M)sA<J»M,M)-SMl 
B < J*M*M) *1 « *2. #SMUl MP 
30  C(J*M»M)sC(J»M*M)“5Pl 

2  CONTINUE 

C.. .APPLY  SYMMETRY  B.C.  IMPLICITLY 
00  4  L= 1 ,4 

fl(2,L»l)*B<2*L*l)*A<2,L*l) 

35  8<2*L*2)s8(2*L*2)*A<2*L*2) 

B(2*L»3)=B(2*L*3)-A(2*L,3) 

8<2,L»4)rB(2,L*4) ♦A(2,L.4) 

4  CONTINUE 

SMlsSMUIMP*0(l«K>/0(2tK) 

40  8(2*1*1)36(2*1*1) -SMI 

8  <2*2*21 *8  <2.2.21 -SMl 
H (2,3*3) XB <2*3*31 *SMl 
8 <2,4.41 cB (2*4.4) -SMl 

C. .. IMPOSE  OUTFLOW  B.C.  USING  LINEAR  EXTRAPOLATION  IMPLICITLY 
45  J=JM 

SP1*SMUIMP*D (J*1.K)/D(J,K) 

00  6  M*1 *4 
DO  S  L=I ,4 

A(J,L*M)=A< J,L«M)-C< J,L*M) 

50  S  H<J,L,M)=8<J,L*M)*2.0*C<J*L*M! 

A<J*M«M) sA ( J.M.M) *SP1 
H<J,M,M)3B<J,M,M)-2.*SPI 

6  CONTINUE 

C...FHL  FORCING  FUNCTION  FROM  S  ARRAY 
55  DO  7  J=?.JM 

00  7  M=1 ,4 

7  EF<J.M)=S<J,K,M) 
return 

F.NU 


1  SUBROUTINE  LBLTRB(J) 

C<>MM0N/C0M1/JMAX*KMAX, JM*KM,XMACH.GAM,GAMM1 *CN*DT*SHU* JCS.PRT* 

1  IPPT  *H, OMEGA, I T  *T  All* 1 TER.tNT .PTORT  *PINF*RINK*0INF,CINF,PT»ITS» 

2  IR1  » IW2, I AFBO* IGEOM, TM* I VIS, II P AN.CE  *CC » JNM.PF Y»PkO»CVIS*CVI SI » 

5  3  T»A*ITWA*LIP»KRFS*SMJIMP,HTINF *FT 1 NF  *S1 NF *E I INF.RE YlN.SUM (40) * 

40ETT (40 1 ,DETL<40) *tT<40) »TH<40) ,TTF,F ACTB.FACTT.REYNLO.PRTURB 
C0MM0N/C0M2/X (40*40) , Y (40*40) * XFX (40* 40*2) * XF Y (40*40 *2) *0(40*40) 
COMMON/C 0M3/0 (40*40*4) *FF (40*4) ,c (40 • 40 *4) ,S (4 I « A8 (4,4) *HVFC<40,4) 
C0MM0N/C0M4/A (40*4,4) *8(40,4,4) ,0(40*4*4) ,HD(40*4*4) * 

10  1U0 (40  *4*4) » AX (40 ) * AY ( 40 ) ,BX(40) ,RY(40) 

DO  1  K*1 *KMAX 

C...LOAO  BLOCK  B  MATRIX  EVALUATED  AT  N  TH  LEVEL  FOP  ALL  K  INTO  HD  APRAY 
CALL  ABMATX ( J»K ,2) 

DO  1  M*1  *4 
15  00  1  L*1 ,4 

1  HO (K*L*M)*AB(L*M) 

C...FII L  OFF-OIAGONAL  AND  DIAGONAL  ELFmEnTS  BASED  ON  A  2-NO  ORDER 

C.. .CENTRAL  DIFFERENCE 
DO  2  K>2*KM 

20  SM1«SMUIMP*0(J*K-1)/0(J*K) 

SP1»SMUIMP*0(J*K*1)/0(J*K) 
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DU  2  M=1 »A 
00  3  L= 1 »A 

A(K,L.M)a-HO(K-l,L»M)*H 

B(K,L»M)aO.O 

3  C«K,U»H)sHD(K*1,L*M)*H 

,  .SET  B  ON  THE  DIAGONAL  REPRESENTING  ThE  IDENTITY  MATRIX  TO  ONE 
A(K.M.M)xA (K.M.m)-SMI 
B(K,M»M)=1.»2.*SMUIMP 
C(K,M,M)aC<K,M,M)-SPl 

2  CONTINUE 

>  •  ADD  SOURCE  TERM  IMPLICITLY.  UD  REPRESENTS  THE  DH/DQ  OF  SOURCE  TERM 
IF ( JCS»£Q.0)GOTO5 
DO  4  K=2 »  KMAX 
DO  4  M= 1  *  4 
DO  4  L= 1 *4 

4  B(K,L.M)=B(K.L.M>*UD(<.L.M) 

5  CONTINUE 

■  .ADD  VISCOUS  TERMS  IMPLICITLY 
IF(IVIS.GT.O)  CALL  VSMATB(J) 

..APPLY  BODY  B.C.  IMPLICITLY  FOR  NOSI  IP  VISCOUS  FLOW 

..FILL  FORCING  FUNCTION  I  ROM  5  ARRAY 
DO  7  K=2,KMAX 
DO  7  M=1 ,4 

7  EF(K.M)=S(J»K*M) 

RETURN 

END 


SUBROUTINE  LUDtC ( A ) 

DIMENSION  A (4.4) 

COMMON/LUD/  Lll.L2i.L22.L31 .L32.I  33.L 41 .L4?,l  43.L44, VI . V2. V3. V4. 
1  U12.U13.U14.U23.U24.U34 

REAL  L11.L21.L22.L31.L32.L33.L41 .L42.L43.L44 

SUBROUTINE  COMPUTES  L-U  DECOMPOSITION  ELEMENTS 
L 1 1  a  A(l.l) 

VI  a  1./L11 
U12  a  V1*A ( 1 ,2) 

Ul  3  a  VI  *A  ( 1 . 3) 

111*,  a  V1*A(1.4) 

L21  a  A (2,1 ) 

L22  a  A (2,2)  -  L21*U12 
V2  a  1./L22 

M2 3  a  (  A (2. 3)  -L21 *U13) •  V2 
M24  a  (  A (2,4)  -  L2 1  * J14 ) *  V2 
L31  =  A  <  3 , 1 ) 

L 32  a  A (3,2)  -  L31*Ul2 
L33  a  A (3,3)  -  L31*U13  -L32*u23 
V3  a  1./L33 

U34  3  (  A (3,4)  -  L31*U14  -  L32*M?4)*  V3 
L41  a  A (4, 1 ) 

142  a  A (4,2)  -  L4 1 *Ul 2 
L43  a  A  (4,3)  -  LA  1*1)13  -  L42*U23 
L44  a  A  (4,4)  -  LAI *Ul  4  -LA2*U2a  -LA3*J3A 
V4  a  l./LAA 
RETURN 
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SUBROUTINE  OUTPUT (L) 

COMMON/COM1 / JMAX  »KMAX ,  JM.KM.XMACH.GAM.GAMM1  t  C.N  , DT  .  SMU  •  JCS ,  PPT  » 

1  IPRT.H.OMtGA.IT.TAu.ITER.FNT.RTORT.PINF.PINF  .OINF.ClNF.PT. ITS* 

2  iki,iw2.iafbd.igeom,tm.ivis.i (pan.cf.cc. jnm.pfy.prd.cvis.cvisi, 

3  TWA.ITWA.L IP.KRES.SMJIMP.hTINF  .FT INF  »SINF »E I INF » Rt YIN.SUM (40) * 
40ETT (40) »DETL (40) ,ET(40) »Th (40) , I TF .F ACTb.FACTT .RtYNLO.PPTURB 

C0MM0N/C0M2/X (40.40)  t Y(40.40> »XF* (40.40,2) ,  Xt:  Y (40-.40  »?) .0(40.40) 
C0MM0N/C0M3/Q (40 .40 ,4) .EF (40.4) .<:(40.40.4) ,G (4) .Ad (4.4) «HVFC(40.4> 
DIMENSION  BHO (40*40) • SL (40 ) .CON ( 8 ) »CP(40) «RCP2(40) « DRAG (40 ) » 
1LP(40).XSL(S00) .YSL(SOO) 

DIMENSION  DCON(40) .EC ON (40) .TPI40.3) ,P(40»3> 

DIMENSION  R0(40,3) .OHO (40, 3) ,W(40»3) 

DIMENSION  PHI (3) *C(3) ,C2(3) ,CPHT (3) .U0I40.3) » VO (40. 3) 

DATA  FLAG/1./ 

GO  TO  (1.2.3.4.S.M.L 
1  CONTINUE 

OUTPUT  FLOWFIELD  DATA 
IF (FLAG.EO.O.)  GO  TO  118 
READ  (5. 1 19)  (LPdS.l3l.KMAX) 

FLAGsO. 

118  CONTINUE 

119  FORMAT (8011) 

SUM (2) S5QRT (X(2»1)**2*Y(2*1) **2) 

DO  11  J=3»JMAX 

11  SUM<J)=SUM(J-1)*SQRT((X(J.1)-X(J-1,1>)**2*(V(J.1>-Y(J-1.1))**2> 
PMS=0.0 
PERRMX*0.0 
KSL=  1 

SUM ( 1 ) *-SUM (2) 

DO  10  K= 1 . KMAX 
IF(LP(K).EQ.0)  GO  TO  131 
WRITE (6.120)  K 

120  FORMAT (*0*. ‘SECOND  INDEX**, 13./ ) 

IF(K-l)  303,304,303 

303  CONTINUE 
WRITE(6.117) 

117  FORMAT (*  1ST*,4X«*P/PINF*,4X,*R0/RINF**4X,*U/QINF**5X,*V/QINF*.?>X. 
•*S/SINF*,4X.*HT/HTINF*.SX»*MACH*,BX,*CP*.9X.*X*»10X»*Y*.7X. 
#*EI/EIInF*> 

GO  TO  309 

304  WRITE (6, 301 1 

301  FORMAT (*  1ST*.4X»*P/PINF*»4X,*  S  *,4X**U/0INF**5X.*V/0INF*.SX, 
•*S/SINF*,4X«*HT/MTINF*,SX,*R/RI*.BX»*CP*.<>X.*X*»10X»*Y*»7X. 

**e i/ei inf*) 

309  CONTINUE 
131  CONTINUE 

DO  F>6  Jsl.JMAX 

EN*Q ( J»K ,4) *0 ( J.K ) 

RMO (J,K)*0(J,K»1)*D(J,K) 

U=0(J»K»?)/0(J,K.l) 

V.*0(  J*K,3)/0(  J.K,  1 ) 

PA*GAMM1*(EN-RH0( J»K)*0.5*(U*U*V*V> ) 

CPP= (PA-1 . ) / (0.S*GAM*XMACH**2) 

ENTRO»PA/RHO ( J. K) **GAW 

HT*GAM/GAMM 1 *PA/RHO ( J.K ) .0 ,S* (U*l  I* V*V ) 

SS=SORT(GAM*PA/RHO(J,K) ) 

Ul«()/0INF 
Vl *V/OlNF 

htisht/htinf 

ElR*(PA/RHO(J»K) )/(GAMMl*EIIKF) 

PERR*ABS(HT-HTINF)*100.0/HTIMF 
lF(PERR.GT.PERRMX)  PERPMX3RERR 
PMS*RMS*PERR**2 
SL ( J) *SORT (U*U* V*v ) /SS 
IF(LP<K) .EO.O)  GO  TO  66 
IF(K-l)  306.307,306 
306  CONTINUE 

WRITE (6, 121)  J.PA,RrtO( J.K) .Ul.Wl .kNTRO.Hll.Sl (J) »CPP»X( J.K) ,Y( J.K) 
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GO  TO  308 

307  WRITE (6*121 )  J.PA.SUM( J) ,U1 *V1 *f nTRO*HT t *RHO( J.K) • 

1  CPP*X ( J*K ) «V(J«K) *E  I R 

121  F0RMAT(I3*11E11.4) 

308  CONTINUE 
66  CONTINUE 

UO  10  J*3* JMAX 

IF(<SL(J).LE.1.0.AND.SL(J-l).6E.1.0>.0R.(Sl  (J) *GE . 1 .0. ANO.Sl (J-l). 

1  LE.1.0))  GO  TO  12 
60  TO  10 

12  JSL= J 

JSLM»JSL-1 

COEF*( 1 .O-SL  <  JSLN) > / ( SL ( JSL )  -SL  ( JSLM ) > 

XSL (KSL) *X ( JSLM.K )*COEF#(X ( JSL  «K ) -X ( JSLM.K) > 

YSL (KSL) *Y (JSLM* K ) *COEF* ( Y ( JSL*K ) -Y(JSLM*K)) 

KSL*KSL*1 
10  CONTINUE 
WRITE(6*111) 

111  FORMAT (*0*»#  SONIC  LINE  LOCATION**/) 

KSLsKSL-1 
00  122  K*1 *KSL 

122  WRITE<6.110)  XSL(K)«YSL(K) 

110  FORMAT <*  XSL**«E11.4.3X.<»YSL»**F1  1.4) 

RMS=SQRT (RMS/JMAX/KMAX) 

WRITE(6* 107)  PERRMX.RHS 

107  FORMAT («0***  PERCENT  ERROR  In  HT=**F 12.4.3X**  RMS  OF  PERCENT  ER 
2R0R  IN  HT»*«E12.4,/) 

T0bM2*2 . /GAM/XM ACH**2 
DO  61  J=1 « JMAX 
RO=Q ( J* 1 • 1 ) *0 ( J* 1 ) 

E=0 ( J* 1 *4 ) *D ( J* 1 ) 

U*Q(J.1*2)/0(J*1*1> 

V=Q ( J* 1 *3)/0(J*l*l) 

PA-*GAMM1*(E-RO«0.5*(U*»2*W**2)  ) 

CP(J)*T0GM2*(PA-1 .) 

RCP2(J|aY(J,l)**2 

61  CONTINUE 
SUM?sCP  <  2 ) *RCP2  <  2 ) 

IF(JMAX-l)  64*63*62 

62  00  65  J*2 . JMAX 
SUMlaSUM? 

SUM2»SUM2*0.5*(RCP2( J)-RCP2( J-l ) )  * (CP ( J) *CP ( J-l >  > 

Rb=Y (JMAX* 1 ) 

65  OP AG (J-l ) sSUMl/RP**? 

63  ORAG ( JMAX ) »SUM2/PB**2 
WRITE (6* 164)  DRAG(JMAX) 

164  F0KMAT(1X«*PRESSURF.  URAG  =**5X,F1  3.  10) 

64  CONTINUE 
RETURN 

2  CONTINUE 

OUTPUT  E  AND  E  CONSERVATIVE  VARIAHI  ES 
WRITE (6* 103) 

103  FORMAT(*0«*37X.*CONSERVATIVF  VAPTAbLFS*) 

00  7  K«1.KMAX 

WRITE (6 *104)  K 

104  FORMAT (*0***K»*» 12 t //»4X**J**6X**t i*» 10X«*F2*. 10X**E3*« 10X**F4** 

2  10X**F1*»10X»*F?*»10X**F3«*10X,*F4*./) 

00  7  J*1 • JMAX 

CALL  EEC0N(J«K«1) 

00  8  N«1.4 

8  CON (N) «G (N) 

CALL  EEC0N(J«K*2) 

00  9  N*1 .4 
NN«N*4 

9  CON (NN) *G (N) 

7  WR1TE(6«105)  J* (CON (N) *N> 1*8) 

105  FORMAT (I5«8E12.4) 
kLTURN 

STORE  DATA  ON  TAPE?  FDR  RESTART 
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140 


145 


150 


155 


160 


165 


170 


175 


180 


185 


190 


195 


200 


205 


C 


3  CONTINUE 

WRITE!?)  (4X(J,K)  ,J*1,JMAX) ,K*1 tKMAX)  t 
1  ( ( Y ( Jt K ) t  J*1 , JMAX) ,K*1 tKMAX ) t 

1  ((D(J,K)tJ*l t  JMAX )  ,K*  1 tKMAX) • 

1  ( I (XEX( JtKtN) • J*Lt JMAX) »K*'l tKMAX) ,N«) *2) t 

1  ( I (XEY ( JtK  »N) ,J*1 t JMAX), K*1 tKMAX) ,N*1 ,2) t 

1  ( ( (Q( JtKtN) ,J«lt JMAX) tKs) tKMAX) tN*l,A> t 

1  JMAX, KMAX tXMACH, RAM, IT t TAU, f ACTTt (OF TI ( J) t J*1 t JMAX) 

WRITE 46,113) 

113  FORMAT (*0*t ‘SOLUTION  MAS  BEEN  STORED  ON  DISK*) 

RETURN 

STORE  INITIAL  DATA  FOR  AFTERBODY  CAL.  USING  NSWC  INVISCIO  30  COOF 

4  CONTINUE 
REAO(6tlOO)  JWRIT 
WR1TE<6,401)  JWRIT 

A01  FORMAT (*0*t*OATA  AT  J*«,I3,2Xt*  |S  STORED  FOR  AFTERBODY  CAL.*) 

J* JWRIT 
NC=KMAX 
MC=3 

Z7=X ( J, 1 ) 

ATTACK*0. 

ACHsXMACH 

YAW*0« 

GAMMA* 1.4 
R I NF  * l  • 

DlNF*l. 

PHIO*3. 141592 
NGAS*0 
NTEST*0 
RRXsO. 

FN*0  . 

F  Y*0. 
f  A=0. 

MX*0 

MY*0 

mZ*o 

FNZ*0. 

FVZ*0. 

FAZ*0. 


MX?*0 

MYZ*0 


mZZ*0 

OPH*?.*ATANI|.) 
00  35  M*1,MC 
J* JWRIT 

PHI IM)*DPM*(M-1) 
C (M) *Y ( J ,KMAX ) 


C«M)*T« J,RM**» 

CZ<M)*(Y(J*1,KMAX)-V( J-l tKMAX) >/(X( J»1,K«AX)-X 4 J-l tKMAX)) 


CPMI (M)*0. 

DO  31  K*1 tKMAX 
JAaJtl 

IF (X(JAtK) .LT.ZZ)  J* JA 
Jl«J*l 

RD (K  tM) *Y ( JtK) 

DR0(KtM)*Q4JtKtl)*l)  <JtK) 

VQ(K.M)=0. 

EN*0 4 JtK ,4) *D  ( JtK ) 

UA*0 ( JtK ,3) /QCJtKtl) 
VA*0(JtKt2)/QCJ,K,l) 

Rl*Y4JltK) 

Dl*04Jl,Ktl)*D  (JltK) 

E1*0( J1»K»4)*D  (JltK) 

Ul*0(JltKt3)/0(JltKtl> 

Vl*Q(JltKt2)/0(Jl,Ktl) 

RATIO* (7Z-X (JtK) )/<X( JltK) -X ( JtK) ) 
UA*UA*RATI0*(U1-UA) 

VA*VA*RAT 10* ( Vl-V* ) 

EN*EN»RATIO* (El-EN) 

DRO(KtM)*DRO(KtM) *RAT 10* (Dl-ORO (K«M) ) 

RD(KtM)*RD(K,M)tRATIO*(Rl-RD(K,M) ) 


V,- 
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UO(K*M)*UA 

W(K,M)«VA 

R(K,Ml* (GAM-1. )*(EN-DR0(K.M) *0.5* <UA*uA»VA*VA) ) 

IF(M.GT.l)  GO  TO  31 

*RITE<6,101)  K.P (K.M) ,DRO(K«M) ,W(K,M),VO(K.M> .UQ(«.M> ,RO(K«M) 

101  FORMAT <*0*»*K»*. I2«2X,*PR£SS**,F1 0.4.2X ,*DENS**.E10.4.2X« 

1*AX  V£L**.£10.4.2X.*CIKCUM  VEL»*,E10.4,2X,*RAO  VEL»*.E1 0 ,4.2X. 
2*RO«*,E10.4) 

31  CONTINUE 
100  FORMAT  < 15) 

IF(X(J*R1T.KMAX).EQ.ZZ)  GO  TO  34 
C (M)*RO(KMAX»M) 

CZ<M1W<Y<J1 .KMAX) -Y(J. KMAX) 1 / (X ( Jl .KMAX) -X ( J.KMAX) ) 

34  CONTINUE 

IF (M.GT . 1 )  GO  TO  3b 

WRITE (6«302)  MtPHl (M) ,C(M)  tCZ(M) ,CPMI (M> 

302  F ORMAT (*0*»I5*4(?X*E10.4) .//) 

35  CONTINUE 
K  =  1 

WRITE (7)  NC.MC* ATTACK* YAW* ACM .GAMMA *R1NF .DINF .PMIO.K.77, 

1  NGAS,  NTfST*  HMX. 

2  FN.  F Y.  FA.  MX.  MY.  M7.  FNZ.  FY7.  FA7.  MXZ.  MY7.  MZZ. 

3  1PM I  (M)  .  C<<4).  CZ  (M)  »  CRMllW).  M  «  1  .  ML). 

4  1 (RO(N.M) .UO(N.M) .VOIM.H) ,W(N.M) .P(N.M) .DRO(N.M) . 

5  M  *  1  .  MC)  .  N  =  1  .  NCI 
RETURN 

OUTPUT  DETAILED  RESIDUE  INFORMATION 

5  CONTINUE 

IF (KRES.EQ.01  GO  TO  33 
WRITE16.U4) 

DO  32  Ks2.KM.KRES 
WRITE  <6.1 151K 

WHITE <6. 116) (J.S(J.K.l) .S1J.K.2) .S1J.X.3) ,S < J.K .41 . J»2, JM) 

32  CONTINUE 

33  CONTINUE 

114  FORMAT (1H0,T4S.*0ETAILED  RESIDUF  INFORMATION*. /.T45.?8(lH*>> 

115  FORMAT  1 1H  «T60»*StCON3  INDFX  »*,14) 

116  FORMATUH  » I4.4E1S.7.2X.  I4.4E15.7) 

RETURN 

OUTPUT  MEAT  TRANSFER  INFORMATION 

6  CONTINUE 
REYsREYIN 

CMU1J*  <TWA**1 .5)*CVISl/(TWA*CVlS) 

TSTAG*1 . .0.5*GAMMl*XMACH**2 
CC«CMUU/ (REY*PRD* 1 TST  AG-TWA) 1 
WRITE16.221) 

221  FOKMAT (*0*.*OISTRIbUT10N  OF  STANTON  NUMBER*) 

WRITE <6.2201 

220  FORMAT <*0*.2X»*J*.15X»«S*»18X»*5T**18X»*T I*. 18X»*T2*.18X.*T3*) 

DO  65  J* 1 . UMAX 

ECON  ( J)  sSQRT  (XEX  ( J.  1  »2)**2.XEV  (.1.  1  .2)  **2) 

DCON ( J) slXEXlJ.1.11*  XEX(J.I «2>*XtY(J.l.i)*XFY(J.l«2) 1/ECON1J) 

DO  69  K*1 .3 

P (J.K 1*10 (J.K, 4 )-0. 5* (Q<  J.K  .2)  •*?♦<! (J»K.3)**2)/Q(J.K.1)1 
1  *0 (J.K 1 *GAMMl 

TR( J.K)*P(J.K)/D(J.K)/0( J.K.l) 

69  CONTINUE 

00  79  J»2. JM 

TNN»DCON ( J) *0 .5* (TR(J»l*l)-TP(J-l»l) 1 .ECON ( J) • (-3.*TP< J* 1 >*4.*TP(J 
1  *2) -TP (J.3) 1*0.5 
ST»CC*TNN 

WRITE <6. 216)  J.SUM(J) tST.TP(J.l) .TPIJ.2) .TP (J.3) 

218  FORMAT ( I5.5E20.6) 

79  CONTINUE 

IF(ITF.EO.I)  PETURN 

DO  89  J»2 . JM 

wRI TE  <6,50 1 1  J , SUM ( J I 

WRITE (6.5021 

DO  89  K* 1 .KMAX 


-• 
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PJK»(Q(J,K.4)-0.5*(Q(J*lC.2)*«2*n(J.*.3>**2>/0(J»K»l)) 

1  *0(J*K)«GAMMl 

200  TJK*PJK/D(J.K)/0( J*K*1> 

SJK*(X(J»K)-X(J*1>  )<*<*2*(Y(J»k)-Y(J»l>  >**2 

SJKsSQRT(SJK) 

l)JKsQ(  J«K«2)/0(  J*K«1 ) 

VJK*Q(J*K*3)/Q(J«K*1 ) 

285  SUV*PJK<*GAM/  (Q  ( J*K  *  1 )  *1)  ( J«K )  ) 

UMACH*SORT  (  (UJK*«*2*V  J<**2)  /SUV) 

WHITE  (F>*503)K  »SJK *  T JX*RJK .UMACh.UJK  *  V JK 
89  CONTINUE 

501  FORMAT  (5X  **«)*•*  IS* 10A»*S/L**»F8,a) 

290  502  FORMAT (8X»4K**SX»#N0H*4AL  UISTANCF**10X,*lFMPFRATURt*»10X. 

1  *H»ESSURE*«10X**MACH  NO.*.10X.*u-vEL0CITY*»10X**V-VtLnCtTY*) 

503  F0kMAT(5X*I3*SX*ClS.b*bX*t:iS.5*5X*klS.5*bX*nS.5«5A*ElS.5*5XtE15.5 
1  > 

RETURN 

295  END 


1  SUBROUTINE  RESIOU 

COMMON/COMl/JMAX*K*AX« JM**M.XMArw*GAM.GAMM) .(  N • OT t SMU* JC«. .P^T  « 

1  IRRT#M.OMFGA.IT«TAU«lTtP»FNTtHTOHTtPlNF.HINF.OlNr *CINF*PT* ITS* 

2  lRl*IW2*IAFBD.IGtOM,TM, I  VIS. I TPAN.CF *CC* JN4.RF  Y.PKO.CVIS.CVISI , 

5  3  TWA* ITWA.LIP.KRfcS.SM JIMP. HT INF , FT  INF tSlNF  «FI  INF • RE YIN* Sum (AO) • 

4DETT (40) *DETL (40) «ET (40) «TH(40) «  )TF *F  ACTU.F ACTT *k£ YNLO.PPTUkR 
COMM0N/C0M2/X (40*40) »Y (40*40) *XFX  (4U*  40*2  )*  Xt- Y  (40*40*2) *0(40*40) 
C0MM0N/C0M3/0 (40 *40*4) *EF (40.4) (40*40*4), G (4) *Ad (4*4) *hVEC (40.4) 
RSOMAX*0*0 
10  RSUTOT*0.0 

01234*0.0 
DO  100  J*2* JM 
DO  100  K*2.KM 
RSDSQR-0.0 
IS  DO  5  L*1 *4 

QLMNT*S ( J  *K . L ) 

RSDSOR-RSDSOR *QLMNT 
IF (OLMNT.LT. 01234)00105 
01234 *OLMNT 

20  J1234*J 

K12344K 
L12344L 
5  CONTINUE 

IF (RSOSQR.LT. RSDMAX)G9T010 
25  RSOMAX*RSOSQR 

JRESDU*J 
KRE S0U*K 
10  CONTINUE 

RSDTOT*RSOTOT*RSDSOR 
30  100  CONTINUE 

RSOMAX*SORT (RSOMAX) 

RSOTOT*SORT (HSOTOT) 

01234*S0RT (01234) 

PERCNT*RS0MAX/RSDT0T#1 00 . 

35  WRITE (6*200) JRE SOU* KRESOU»RSDMAK,RSOTOT .kERCNT* J1 234 .K1234.L 1234* 

1  01234 

200  FORMAT (*  RESIDUE  INFORMATION»»9X.2I5.3F10.5**S(**I3«4*4*I3*4»**T2, 

»  *)**»F10.5) 

RETURN 

40  ENO 
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1  SUBROUTINE  WHS 

COMMON/COM1/JMAX.KMAX. JM.XM.XMACM.GAM.GAMMl.CN.DT.SMU. JCS.PRT. 

1  IPRT.H.OMFGA.IT.TAu.ITFR.FNT.RTORT.RlNF.klNF.OINF.ClNF.kT.nS. 

2  IR1  .  IW2.  I  AFBO*  IfiF.OM*  TM.IVIS.I  iPAN.Ct  .CC«JNM»RFY»PkD»CvIS«CVISl» 

S  3  TWA. ITWA.LIP.KRES.SMJIMP.hTINF «FTINF*SlNF»ElINF»RtYIN»SUM(40)» 

4DET  T  (AO )  »DETL  (40)  .t T  ( 40)  »TH  (40 >  ,  TTF  *F  ACTb.FACTT.RtYNLD.PRTURB 
C0HM0N/C0H2/X (40»40) .Y (40.40) .XF*(40, 40 .2) , XFY (40.40.2) .0(40.40) 
C0MM0N/C0M3/Q (40.40 .4) .FF(40.4) .<M40.40.<») .6(4) .Aa(4.4) «HVFC(40*4) 
C...  DATA  Cl .C2.C3/1.0.-1. 0.0.0/  FOk  ?  POINT  ONE SI Of D  DIFFERENCING 
10  C...  DATA  C1.C2.C3/1 .5.-2.0.*0.5/FOk  3  POINT  ONE S I DF D  DIFFERENCING 

DATA  Cl .C2.C3/1 .0.-1 <0.0*0/ 

C...THIS  SUBROUTINE  COMPUTES  THE  RIGHT  HAND  Slot  OF  THE  DELTA  FORM 
C... EQUATION 

C...FOR M  £  CONSERVATIVE  VARIABLES  AND  DIFFERF NCF.  STORE  IN  ThE  S  ARRAY 
IS  UO  1  K»2,KM 

DO  ?  J*1 * JMAX 
CALL  EFCON(J.K.l) 

DO  2  N*1 .4 
2  EF(j,N)=G(N) 

20  C.. .CENTRAL  DIFFERENCE  E  CONSERVATIVE  VARIABLE 

DO  1  N*1 .4 
DO  1  J*2«JM 

1  S(J,K.N)s(EF(J*l.N)-tF( J-l.N) )*H 
C...FOWM  F  CONSERVATIVE  VARIABLES  ANU  DIFFERENCE.  ADO  TO  PREVIOUS  S 
25  C... ARRAY 

DO  IS  J*2.JM 
DO  4  K*1 »KMAX 
CALL  EFCON( J.K.2) 

DO  4  N*1 *4 

30  4  EF(k,N)*G(N> 

C... CENTRAL  DIFFERENCE  F  CONSERVATIVE  VARIABLE 
00  3  N*1 ,4 
DO  3  K*2.KM 

S( J,K.N)«-S(J.K,N)-(EF(K*l.N)-fcF(K-l,N) )*H-HVEC(K.N) 

35  3  CONTINUE 

15  CONTINUE 

C.. .COMPUTE  TURBULENT  VISCOSITY  COEFFICIENT  IF  NFCF SSARY 
C  IF (IVIS.EO. 1 .AND. ITURB.FQ. 1 )  CAM  MUTUk 
C...A0D  VISCOUS  TERMS  TO  RIGHT  HAND  SIDE 
40  IF ( IVIS.EO. 1)  CALL  VSRHSH 

RETURN 
ENU 


1 


S 


10 


15 


20 


SUBROUTINE  SHAPfc 

COMMON/COM1/JMAX.KHAX. jm.km.xmach.oam.gammi .CN.DT.SMU. JCS.PRT. 

1  1PRT. H, OMF GA. I T. TAU. I  TER. t NT. RTOkT.RlMF. PI NF. OINF.CINF.PT. ITS. 

2  Ikl , IW2. I AFBD. IGEOM. Tp* I VI S. I 1RAN.CE .CC. JNM.RE Y.PHO.CVIS.CVIS1 . 

3  TWA.ITWA.l IP.KPf S.SMJIMP.HTINF .FTINE.SINF .El INF.«tYlN.SUM(40) • 
40ETT (40) .OF  TL (40) .£ T (40) .Th(40) , ITF .F ACTH.FAfTT.REYNLO.PRTURM 

COMMON/HOTH/XOl, X02.X03.X04.YDI.Y02.V03.Y04.SL1.SL2.SL3. Pi. R2. 
lk3.R4.CTl  .CT2.CT3.CT4.CT5.CT6.XOO.RBODY 
COMMON/XYPS/Xl .X2.X3.X4.X5.XG.A7.V1.V2.Y3.Y4.Y5.YG.Y7 
C  THIS  SUBROUTINE  READS  AND  wRITtS  CONTROL  PARAMETERS  FOR  NOSETIP  S 
C 

RE AO (S. 1 21 )  Xl.X2.X3.X4.X5.X6.X7 

121  FORMAT (BE  10*0) 

C 

122  FORMAT (20X.10F10.5) 

C 

READ (5.121)  Y1.V2.Y3.Y4.Y5.Y6.Y7 
C 

WRITE (6*131) 

131  FORMAT (•()•»  5X.*X  AND  V  VALUFS  FOR  THE  CONTROL  POINTS*) 

WRITE (6. 122)  X1.X2.X3.X4.X5.X6.X7 
WHITE (6. 122)  Yl.Y2.Y3.V4.Yb.V6.Y7 
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READ (5, 121 )  R1.R2.R3.R4 
WRITE (6*132) 

132  FORMAT(*0*  »5X»*PA0IUS  FOR  CIRCiRAR  ARCS*) 
WRITE (6,122)  R1»R2»R3»R4 

READ (S* 121 )  SL1.SL2.SL3 

WRITE <6, 133) 

133  FORMAT (*0*.5X.*ANGLES  FOR  STRAIGHT  LINES*) 
WR1T£(6,122)  SL1.SL2.SL3 

l)TR«3. 14159265/180. 

Y01>CU 

XOlaRl 

SL1*0TR*SL1 

SL2*0TR*SL2 

SL3«DTR*SL3 

X02aX2*R2*C0S(SLl) 

Y02*Y2-«?*SIN(SL1) 

X03=X3-R3*C0S(SLl) 

Y03*Y3*R3*SIN(SL1 ) 

X04»X5*P4*C0S (SL2> 

Y04sYS-R4*SIN  (SL2) 

WRITE  16* 134) 

134  FORMAT (*0*»5X,*CENTERS  FOR  THE  CIRCULAR  ARCS*) 
WRITE <6. 122)  X01 »Y01 »X02»Y02*X03»Y03»X04»Y04 
RC0NE*Y6/C0S (SL3) 

XOG*X6.RCONE*SIN(SL3> 

OMEGA*XOO 
wRITE(6,422>  XOO 

422  FOHMAT  ( *0*»*XOOaNtW  OWE6A«*.F10.4> 

PH00Y«X00 

CT1  *ATAN(Y1/(X00-X1) ) 

CT2  *ATAN(Y?/(X00-X2)) 

CT3  *ATAN(Y3/(XOO-X3)) 

CT  4  *ATAN(Y4/|XOO-X4>> 

CTS  *ATAN(YS/(X0U-X3)) 

CT6  «ATAN( Y6/ (XOO-Xb) ) 

WRITE (6.135) 

136  FORMAT <*0*»*THETA  VALUES  FOR  CONTROL  POINTS*) 
WRITE (6,1 22)  CT1 ,CT2,CT3,CT4,CIS,CT6 
RE  TURN 
ENO 


SUBROUTINE  SHOCK 

COMmON/COMI/JMAX.KRAX, JM,KM,XRArM,bAM,OAMMl,(  N*0T.SRU, JCS.RBT* 

1  IRRT.M,OMEGA,IT.IAt), ITER,ENT.rTOkT»R1MF,rINF.OInF»CINF,RT,ITS» 

2  IR| » IW2» 1  AFHf), lot  OR. IM, IVIS, IIPAN.CF.CC.UNR.RFYfRRD.CVTS.CVlSl, 

3  TWA, ITVA.L IR»KRE  S.SMJlMp.HTINF ,rT INF.SINF.EI INF, RE VlN, Sum (40) * 
4DETT  (40)  ,l)t  TL  (40)  »fcT(40)»TH(40)  » 1 TF  »F  ACTt»»F  ACTT ,R£YNLO»RRTllRB 

COMMON/COM2/X (40*40) »Y (40,40) ,XFX (40,40,1) »XFY (40*40, 2) *0(40,40) 
C0MM0N/C0M3/0 (40,40,4) ,FF(40,4) ,C (40,40,4) ,0(4) .A* (4,4) ,HVFC(40,4) 
DIMENSION  R(40,3) *RXl (40) ,RFTA (40) .0(40, J) ,0X1(40) *I)E  TA (40) , 

0  W(40,3) ,VXI (40) . VE  TA (40) ,R(40,3) ,RTAJ(40) ,DTS(40) *XST(40) ,YST(40) 
DATA  XST,YST/40*0. 0,40*0.0/ 

C. ..COMPUTE  THE  FLOW  VAR1 A6LES  OnE  ME  Sm  INTERVAL  BFLOW  SHOCK 
RMSaO.O 
OSEm«0.0 
JMM.JMAX-2 
KMMsKMAX-2 
00  3  K*1 ,3 
00  3  J*1 , JMAX 
KK«KMAX-3*K 
7*1 .0/0(3, KK, 1 ) 

R(J, K)*0(3,KK,1)*D(3«AK) 

U(J,K)*0(J,KK,2>*2 

V(J»K)*b(J*KK,3)*Z 
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E2*0(J.KK.4)*D( J.KK) 

3  P(J,K)»(E2-0.5*H( J*K)*(U(J»K)**?*V( J.  *>**?) )*GAMMl 
..COMPUTE  P-XI.  U-Xl.  P-ETA*  U-ETA,  AND  V-ETA  DERIVATIVES 

DO  4  J«?.JM 

PX1(J)»(P(J*1.3)-P(J-1.3))*0.5 
UXI (J)»(U(J*1.3)-U(J-I,3))*0.S 

4  VXI(J>«(V(J*1.3)-V<J-1.3))*0.5 
PXI (i)a-PXI (2) 

UXI(1)«-UXI(2) 

VXI (1 )*VXl (2) 

PXI(JMAX)«(3.0*P(  JMAX*3)-4.0<»P(  IM.3) »P ( JMM, 3) ) *0.5 
UXI (JMAX)«(3.0*U(JMAX,3)-4.0*U< JMf3)*U(JMM.3»)*0.5 
VXI ( JMAX) * (3.0*V( JMAX»3) -4.0*V ( JM*  3) * V ( JMM, 3F) *0.5 
DO  5  J*1 * JHAX 

PETA  (  J)  ■  (3.(I*P  (  J«3)  “4'.0*P<  J»?)  *P(J»l))*0a5 
UETA(J)s(3.0*(l(J«3)-4.0*U(Jt?)*ll(Jal))*0aS 
VETA ( J) * (3,0*V( Jt 3) -4.0*V ( J,2) *V( J*1 > ) *0.5 

5  CONTINUE 

lE(IT.EO.ITER)  vR1TE(6,100* 

DO  10  J* 1 • JMAX 
KsKmAX 
XET*0a 

|)MARaXET*U(J*3)«XEA(J,K.l)*V(J,3)*XEY(J»K,l) 
VHAR*XET*U(J»3)*XEX(J,K»2)*V(J,3)*XEV(J,K»2) 

RCS*GAM*P ( J*3) 

.. DETERMINE  SHOCK  TIME  STEP 

SPSnD«SORT ( 6AM*P ( J, 3) /R ( J,3) ) 
f  TAT*-(XEX ( J,K,2) *XST ( J) «XEY ( J«K*2) *YST ( J) > 

SlGA*ABS (IJtlAR) *SPSNO*SORT (XEX ( J,K  *  1 ) **?*XFY ( JtK *  1 >  **2) 

SIGR*ABS (ETAT*VBAR) *SPSND*SQRT ( xFX ( J»<,2) **2*XFY( J»K»2) •*?) 
SlOARsAMAXl (SIGA.SIGH) 

DTS(J)*.90/SIGAts 

IF  (ITaEOalTFR)  vRITt  (6a  105)  J,STOA,SI5P,UTS(J) 

105  FORMAT («  J**,I2«3X,*SIGA>*,E13.5,JX«*SIOH**»F 13.5,JX,*DTSaMUE13.5> 
H*“RCS* (UX I ( J) *XE* (J»K»I)*VXI(J)*XEV(J»K»1)*UFTA ( J) *XEX (  JHt^ZTv, 
l  ♦VFTA(Jl*XEY<J,K,«!)*(Y<J»3>/t(  I.KH) 

DETERMINE  PRESSURE  *T  SHOCK  EXPLICITLY 
1 1  PETA  ( J) *P I J.3» ♦OTS( J) • I-UbAR*Px 1 1 J) -V3AR*PETa ( J)  > 

10  CONTINUE 

C. .. Fill.  HOUNOARY  POINTS  FOR  PRtSSIlRE 
PETAIDaPETAI?) 

PETA(JMAX)*?.0*PE TA(JM)-PETA!JM-1 > 

C.. .SMOOTH  PRESSURES  AT  SHOCK  USING  FOURTH  ORDER  SMOOTHING 
SMUS=>0.S 
DO  14  J*3» JMM 

14  PXI  (J)spfTA(J)-SMUS*0.12S*(PETA(.)-2)>4.0*PETA(J-l>-»6.0*PFTA(J)<»4.0 
>  #PFTA (J*1)*PFTA(J*2)) 

PXI (2)*PETA (2)-SMUS*0. 125* (2.0*PFTA (2)-3.0*PFTA (3) *PETA (4) > 

PXI ( 1 ) >PX I  (?) 

PXI(JM)mpETA( JM)-SMUS*0.12S*(P£TA(JM-2)-4.0*PETA(  JMM) *5.0*PETA( JM) 
»  ~2.0*PETA (JMAX) ) 

PXI (JMAX )S2,0*PXI I JM) -PX 1 1 JMM) 

DO  1  J»1 ,  JMAX 

..DETERMINE  SHOCK  angle  OELTA»ARCTAm(-ETAY/FTAX> 
IjELTA«ATAN(-XEY(J«K,2)/XFX(J,K,?) » 

S0«SIN (DELTA) 

CD»COS (DELTA) 

UlT*QINF*CO 
P2«PXI (J) 

IF (P2«LF .0*0)  GO  TO  6 
Z«GAM*1.0 

XMX*SQRT (0.5/GAM* (P2/PINF*Z*GAMMl ) ) 

QS«CINF*XMX-U1T 
PH»P(J»3) 

RB*R( J»3) 

UB-II(J,1) 

VH«V(J»3) 

EB>PB/GAMM] *0.5*RH* (Ud**2* VB**2> 

U2T  >2.0* ( 1 aO-XMX**2) •CINE / ( (GAM* | .0) *XMX) *U1 T 


L 
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;S 

>: 


R2«R1NF* (P2/PINF*GAMMl/7 ) / ( 1 ,0*GAMMl/?*P2/PlNF) 
U2»0INF«SD**2*U2T*CD 
95  V2*OINf*SD*CO-U2T*SO 

E2=P2/GAMM1*0.S«P2*(U2»»2*V2**2) 

C... COMPUTE  PTAU 

PT Au  C  J) *  (P2-PB) /OTS  < J) 

C... COMPUTE  CONSERVATIVE  VARIABLES  AT  SHOCK 
100  K*KMAX 

DI*1.0/0(J.K) 

G( J»K»1 )*R2*DI 
G(J.K,2)»R2*U2*DI 
0(J«K,3)*R2*V2*0I 
10S  Q ( J,K»4) *E2*0I 

C. . aOETE hhINE  ANGLE  OF  Xl»C0NST  LINE  WTTH  X-AXIS 
K*KMAX 

IF < AHS(X£Y(J»K.l) )-0. 00000 1)  7,7.6 
7  THETA*1. 57079633 
110  GO  TO  9 

6  CONTINUE 

THF.Th*ATAN(XEX(J»K»1)/XFY(J»K»1)) 

9  CONTINUE 

C... COMPUTE  SHOCK  SPEED  IN  X  AND  Y  DlRFCTIONS 
11S  HETA*THE  TA-Ufcl TA 

»jSE=OS/COS(HETA) 

IF  (AHS(QSE)  .ge.  AhSOSFM) )  JOS=J 
IF(ABS(OSE)  .Gt.  AHSOSEM)  )OSEM=OSt 
RMS*RMS*0SE**2 

120  AST ( J) *-OSE*COS (THtTA) 

YST(J)sQSE*SlN( THETa) 

THtTA»THFTA*57.29S /« 

OELTA*OELTA*S7. 29576 
HE"  TA*BETA*57.2957H 

125  IF ( IT .EG, ITER)  WPITEI6.101)  J.ThfTA, DELTA, BETA, XMX,U1T.U?T.0SE. 

>  XST(J) .YST(J).PH.P2»RB.R2»UB,U?.VB»V2.Fb,E2 

>  .PTAU(J) 

C... PROPAGATE  SHOCK 

X(J,K)«X(J,K)«XST(J)*3T 
130  Y(J,K)*Y(J»K)*YST(J)*DT 

C... ADJUST  OTHER  GRIO  POINTS 
XH«X(J,1> 

YB*Y ( J, 1 1 
DXX*X(J.KMAX)-XB 
135  OYYaY  ( J.KMAX) — Yb 

DO  2  K«2,KM 
ETAaET(K) 

X (J.K)  *  XH  ♦  OXX*ETA 
Y(J.K>  »  YB  ♦  OYV*ETA 
1*0  2  CONTINUE 

1  CONTINUE 

BMS=SORT (RMS/ JMAX ) 
mRITE(6. 102)  RMS, JGS.OSEM 
100  FORMAT (*0*»*FR0M  SUB.  SMOCK*) 

145  101  FORMAT (*0*»*J**» I2»4X»*THE TAa*»eI 0.4, lX»*OELTAa*,£l0,4»lX»*BETA 

>  E10.4,/,9X.*MX>*,E10.4,4X,*UllB*«E10.4,3X,*U?T**,tl0.4*?X«*OSE 
•  E 10.4, 2X.*XST»*»E 10.4, 2X »*YST«*»E10,4»/»9X*1 IF 10 .4) 

102  FORMAT (•  RMS  Of  SHOCK  SPEED**, El ?.4»3X.*J**. I3,3X,*MAX  SHK  SPO* 
»  E12.4) 

150  RETURN 

6  CONTINUE 
K«KMAX 

WRI TE (6, 1 03)  J.P2,P(J,3) ,PT AU( J) 

WRITE (6,104)  UBAR, VHAR.PX I (J) ,UX T ( J) , VXl ( J) ,PFTA ( J) ,UFTA( J) « 

155  >  VETAI J) ,RCS,XEX ( J.K, 1 ) , XEX <J,K»?)»XEV(J»K»1) , XFY (J«K,?),V(J,3) 

>  Y (J.K) 

104  FORMAT (5E 15.5) 

CALL  OUTPUT (1) 

CALL  EXIT 

160  103  FORMAT (•  NEGATIVE  PRESS.  AT  SHOCK,  J**» I2.3X,*PN»*»E10.4,3X, 

>  *PO»*,E10.4,3X,*PTAU«*.E10.4) 
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SUBROUTINE  Th>lH  (A*h*C*X*E  * NL *  Nn ) 
DIMENSION  A(2).H(2).C(?).X(2).E  C>> 
X  (Ml  )  =  C(UL)/H(NL) 

F  (NL )  *  F(NL)/H(NL) 

ML  P 1  *  NL  ♦  1 

DO  I  J  :  NLP1*  NU 

2=1./  (b( J)  -A(J)«X(J-1) > 

x(J)  x  C(J)  *2 

f(J>  x  (E*(J)-A(J)*KJ-1)  )*Z 

NUPUL  =  NU  ♦  NL 

MO  2  J1  =  NLPl *  NU 

J  =  NUPNL  -  J1 

r  (J)*HJ)-X(J)«F(J*I) 

RETURN 

END 


SUBROUTINE  VSMATH  ( J) 

(  UMMO.N/COMl/JMAX.NMAX*  JM,K.M,XMArH*GAM*GAHMl  *f  N.OT  *SMU* JCS»PRT « 

1  IPPT  »H»OM.E  GA»  IT*  TA||«  ITFP.t  NT  .PTOKT  »P1NE  *RINE  *OINF«CINF«PT*ITS* 

2  I  Ml  .IW/.IAFhO.IGEOM.TM.IVIS,  I  IRAN.CX  .CC  «  JNM.  kF  Y.  PkD*CVl S  »CW  I  SI  . 

3  TxA.ITWA.LIP.KPFS.SMJlMP.hTINE  ,FTINE  *SINF«E1  iNF.Rfc  YIN*S|)M(40>  • 
4DtTT<40) .OETL (40> .1 T (40) ,TH (40) , f TF »f ACTh.FACTT.RtYNLn.PPTURB 

COMMON/COM2/X (40*40) *  Y ( 40 • 40 ) *XFX (40*40*2) *XE Y (40*40*2) *0(40*40) 
COMMON/COM3/Q  (40*40*4)  ,FF  (40 .4)  ,«  (40 * 40 .4 )  .5  (4)  «  A3  ( 4  *4 )  *HVFC(40,4> 
COMMON/COM4/A (40*4.4) *H (40*4*4) *P (40.4*4) «H3 (40*4*4) • 
lliU  (40  *4. 4)  *AX  (40)  *AY(40)  ,BX  (40)  *t>Y  (40) 

COMmON/VISC/U(40) *V (40) *C1 (40) *C? (40 ) »C3 (40) .04(40) «CS(40) ,CM40) * 
107(40)  *TC  (40)  ,CSI  (40)  *CS2(40>  .LS3(4U>  *CS4(40)  .CSS  (40)  .CSX,  (40)  * 

2CS7 (40) ,PR(40> 

COMMON/ V I SK /CHUK AP ( 40 ) • TUMMU (40.40) 

UATA  PPTk*FPT.03*T3/l.,l .333333373333*. 333333333333*. f>x,6M»(M>6X>6A,6/ 
C...SET  UP  CONSTANTS  NEtUtU  FOR  ADD  INC  VISCOUS  TERMS  OK  S  AND 
C...T  MATpJCfS  IMPLICITLY 
HkE=O.S*DT/kEYNLl> 

GPk=GAM/PRO 

C 

00  10  K=l *KMAX 

C...AOO  non-axisvmmetmic  viscous  terms  of  s  mai pi x  implicitly 
C  these  terms  AHE  OF  SECOND  DEklVAITVt  TvPt 
Pl=1.0/O(J*K,l> 

U(M=0<J.K*?)*R1 

V(K)=Q(J,K,3)*ki 

TT=(0(J,K*4)*P1-O.S*( J(K)*«2*V(k>**2) )*GAMMl 
CMUKAP(K)  =  (TT**1.S)»CVIS1/(TT*C>/TS) 

(iMU=CMIIKAP(K)  *TUMNU(  J»K> 

GKAP=CMliKAP  <K>  *TUHMU( J*K) *PRTUPM 

(,PMKx(,Pk*GKAP 

OY=l./Y(J,lO 

DJaCxmRF /0( J*K) 

GMU.)AC=GPU*()JAC 
FY=XEY ( J.K.2) 

E'X  =  XEX  ( J«K*2) 
f  YS=EY*F Y 
FXS=EX*FX 
EX  YxEX*EY 

Cl  (K)xGm()JAC*(F«T*EXS*E  YS) 

C2 (K) =GMUJAC*EXY*03 
C3(K)  xGf'UJAC*  (EXS*E  RT»E  YS) 

C4(K)xGPRK*UJAC*(tXS*EYS) 

PR(K)xRi 

TC (K) «Q( J»K  *4) *R1- (ll(K) **2*V(K)**2) 

CS1 (K)x&MUJAC*OY 
Cl  ( K ) *C 1 (K)*RR(K)*2. 

C? (K) »C2 (F ) *RP (X) *2. 

C3(K)xC3(K)*RP(K)*2. 

C4<K)«C4<K>*RR(K>*2. 

CS (A ) *T3#CS1 (K)*RR(K)/Y(J*K> 
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ss 


60 


65 


70 


75 


00 


85 


90 


95 


100 


105 


110 


115 


120 


C6  oo  «ev*ojac*rp ( k ) /r < j»k > 

C7  IK) »FX*0JAC*KH (K ) /Y  < J*K ) 

CS2 (K)  »C5  (K)*V(K)*tX 
CS3<K)=C5<K)*V(K)«tY 
CS4<K)=CS<K>*<U(K)«E  X*V(K)*E Y> 

CS5 (K) *C5  <K ) *EX 
CS6 (K )  »C5  (M ) *F Y 

CS  /  (K)  =-GPPK*DJAC*XE Y  (J*K,2)*WK(K)/Y  <  J*K ) 

10  CONTINUE 
C 

OO  ?0  K=2» KMAX 
KP*K-1 

Hl)<K,2*l)=-<Cl  (KP)*U<<P)*C2(KP)*V(Rk>  ) 
hO(k,2,2)=C1  <K«) 

H0<K.2«3)=C?(KH) 

Hi)  (K  «  2*4 ) -0 • 

H0<K*3*l)=-(C2<KP)«U<<P)*C3<Kf<)*VtKR) ) 

HO(K,3*2)=C2<KH) 

Hi)  (H  *3*3)  *C3  CKR) 

Hl>  (K  *3*4)  =0* 

HI)  (K«4«l)s-(  C4  <KP)*fC(KP)*Cl  (Kk)*U(KH)**2*2.*C?«KH)*U(KR)*V(RP)*C3 
1 <KH)*V(KP)»*2) 

H0<K,4,2)=-<C4(KP)-C1 (MR) >*U(KK>  *C2<KR>*V(KH) 
H0(K,4*3)=-<C4<KP)-C3(KR) >*V(RH> ♦C2(KR)*U(KR> 

Hl)<K,4,4)=C4<KH) 

IF (K.f O.RMAX)  GO  TO  20 
KP=K*1 

nU(K.2*1)=-<C1 (KP)*U<<P)*C2tKP)*V<KP) ) 

UO (K  *2*2) =C1 <KP) 

UlMK,2,3)=C2<KP) 

UD(K«2«4)=0. 

UD<k,3,1)=-<C?<KP)*UUP>*C3<KP)*V<KP) ) 

UD(k,3.2)=C2(KP) 

UU(k*3*3)=C3(KP) 

IlL)  (K  *  3*4)  =0  • 

llO(K,4.1)s-(C4(KPI*TC(KPl*Cl (RP)*U(KP)**2*2.*C?(KP1«U(KP)*V(KP>*C3 
1 <KP)*V«KP)**2) 

UD<K,4,?)s-<C4(KP)-Cl  <KP>  ?*U(KP)*C2(KP)*V(KP) 
ni)<K,4,3)=-(C4<KP)-C3<KP>>*V(KP»*C2<KP)*U(KP> 
tlD(K,4.**)sC4<KP) 

20  CONTINUE 

UO  TO  Ks?,KM 
KR=K-1 
EP=K*1 
OO  31  N=2,4 
DO  31  Ms  1*4 

7 (R,N.H)=A(R*N*M)-Hn<R.N*M) 

H(6,N*H)=tUK,N«M)  ♦  HUfKP*N*M)  ♦|iU(Rk»N*M) 

C (R,N.M)=C(K,N*M)-UO(<*N»H) 

31  CONTINUE 
30  CONTINUE 

C...A0O  ADDITIONAL  AX1SYFME  TPIC  CONTHIPUTION  TO  S  MATRIX  IMPLICITLY 
C... these.  TERMS  APE  OF  FIRST  OFPIVATIVF  TYPE 
DO  40  K  =  l.i\MAX 
)'0  (k  .2*  1 )  =CS2 (K ) 

HO (k ,2*2) =0. 
hO (K ,2*  3) =-CS5 (K ) 

HO (K , 3, 1 ) =CS3 (K ) 
nl)<K»3*?)  =0. 

HlMh  *3*3) =-CS6 ( K ) 

HU(K,4*l)=CS4(6)*2.*V(h) 

Hl>  (M  ,4*2)  =-CS?  (lO 
til)(K*4,3)s-C54lK)-CS3(M 

40  CONTINUE 

DO  4 |  Ks?,KM  . 

DO  41  Ns?  ,4 
DO  41  Ms 1,3 

«  (A  *N*M)  sA  <K  *N*M)  ♦HOU-l  ,N,M) 
f  !(V,N,M)sC(K,N*M)-H|)I<*l,N,M) 

41  CONTINUE 
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C...AOO  THL  TERMS  OF  AX1SYMME TRIC  MATPTA  T  IMPLICITLY 
C... FIRST  Al)L>  TERMS  OPE  RAT  I  SO  ON  0  VECTOR 
00  (SO  K=2.KM 

IPS  KP=K*1 

KR=K-1 
I*V=V  <rcp> -V 
IH|XU(M>)-U(KM) 

hi)  (K  .2.  1 )  =-  (CA  (KP)  *1)  (<P)*C7(KP)*V(KP)  ) 

130  H0(k»2»1)=-(C6(KP)*U(KR)*C7(KR)*V(KW) ) 

Oil  <K  ,2,2) =Co (KP) 

Of)  <K,2.2)=C6(KK) 
hO(K»2»3)=C7(KP) 

UO <K  .2 » 3) =C7 (KR ) 

13S  h0(k.2.4>=0. 

UD(K.2.4) =0. 

H0(K,3.1)=-2.*C6<KP)*V(KP) 

HLMK.3»1)=-2.*C6(KH)*VIKR) 

mO(K.3.?)=0. 

140  UD(K»3»2)=0. 

Hb(K.3»3)=2.*C6(KP) 

HD<K,3»3)=2.*C6(KR) 
lif)  (K , 3,4 ) =0. 

HOIK. 3*41=0. 

14S  m0(K,4,1)=CS7(KP>*TC(<P)-C6<KP)*(U(KP)**2*FRT*V(KP>**2>-C7(KP) 

♦  *b(KP)*V(KP)/3. 

U0<K.4. 1)=CS7(KR)*IC(<P)-C6(KR)* (U (KR) **2*FRT*V (KR) *•?) -C7 (KR) 

♦  *U(KR)«V(KW)/3. 

mU(K,4»2)=-(CS7(KP)-CS(KP> )«U(KP)-T3*C7(KP)*v(KP> 

ISO  Ub(K,4.2)=-(CS7(Kk)-C(>(KR) )»U(KP)-T3*C7(KR)*V(KR) 

hD(K.4.3)=-(CS7(KP)-FRT*CS(KP) 1 *V IKPl *C7 (KP) *U IKP) 
Ub(K,4.3)=-(CS7<KR)-FRT*C6(KR) ) *V (KR) *C7 (KP) *U (KR) 

Ml)  (K  .4.4  )  =CS7  IKP) 

U0(K,4.4)=CS7IKR> 

155  DO  M  N=2.4 

00  A1  M*1.4 

A(K,N.M)=A(K,N,M) *UO(<.N.M) 

M  C (K ,N»M) =C  (K  ,N,M)  -E10  ( <  »N»M) 

H(K,3,ll=H(K,3,l)*h.*C5(K)*V(K) 

1 AO  R(K,3.3)=B(K.3.3)-A,*CS(K) 

RIA.4.1 )=H(K.4. 1) ♦4.*CS(K)*VIK)**2 
«IK.4.2)=dIK.4.2) ♦5.*C7(K)*0V/3. 
H(K,4.3)=b(K»4.3)-5.*C7(K)*0U23.-4.*C5(K)*V(K> 

60  CONTINUE 

IAS  RETURN 

tNi) 


1  SUBROUTINE  VSPHSR 

C0mm()N/C0M1/JMAX.»>MAX,  JM.KM.XMACH.GAM.GAMMl.CN.OT.SMU.ICS.PHT. 

1  IPRT.H.OMtGA.IT.TAU.ITEP.ENT.RTORT.PlNK.RINF .OINF .CINF .PT. ITS. 

2  lRl.Ih2.IAFB0* IOtOM. IM.IVIS.IlPAN.CF ,CC»  JNM.RF'Y.RRD.CVIS.CVISI . 

5  3  TkA.ITWA.LIP.KRES.SMJIMP.HTINF .FT  INF .S INF ,£ I  INF • RE  YIN, SUM (40)  , 

4l)E  TTI40) .0ETLI4  0) >E  T (40) .THI40) ,TTF ,E  ACTh.FACTT.REYNLO.PRTURH 
C0MM0N/C0M2/X  (40.40)  ,Y  (40.40)  .AFX  (40 , 40 .2)  •  XF  Y  (40 , 40  .  2 )  .D  (40 .40 ) 
COMM0N/C0M3/0 (40,40.4) .EF (40,4) ,c (40 .40.4) ,G (4) • Ad (4.4) ,HVFC(40,4) 
C0MM0N/C0M4/A (40.4,4) ,6(40,4,4) .C(40,4,4) ,HD(40,4,4) « 

10  luO (40,4.4) ,AX (40) .ay (40) ,HX (40 ) ,RY (40 ) 

COMMON/V I SC/U (40) »V (40 ) .Cl (40) »C? (40 ) ,C3 (40 ) ,C4 (40 ) ,C5 (40 ) »C6 (40) . 
107(40) »TC(40) ,CS1 <40) ,CS2 (40) »CS3<40) ,CS4 (40 ) ,CS5 (40) »CSA(40) , 

2CS7 (40) , KM ( 40 ) 

CC'MMON/V I SK /CMUKAP  ( 40 )  ,  TUHmII  (40.40) 

15  DATA  PRTR.F RT .C3.T3/1 • • 1 . 333333333333. . 33333333333 J, .6A6AA66AA6A6/ 

DO  1  J= 1 «  JmAX 
00  1  K=1,KMAX 
1  TURMU( J.K)*0. 

HHE*0.5«OT/HEYNLO 
20  tiPR*GAM/PRO 

DO  30 
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00  .31  K=1  .KMAX 

Pl*l  • 0/0 (J»K*  1 )  ' 

RR(K)=R1 

U(K)=Q(J,K,?)*S1 

V(K>=Q(J,Kt3)*Rl 

TT*(Q(J,K.4)*R1-0.S*( J(K)**2*V(k>**2) )*GAMMl 
CMUKAP(K)s(TT**1.5)*CVISl/(TT*Cv|S) 

GMU=CMUKAP(K> ♦TURMU(J.K) 

6K  AP=CMIJK  AP  (K)  ♦  TURMU  ( J,  K )  *PRTUHH 

gprk=gpr*gkap 

0V=1./Y(J*K) 

0JAC=HRE/0< J.K) 

GMUJACsGMU*OJAC 

£Y=XFY(J*K*2) 

EX=XEX(J,K»2) 
t YS=EY*EY 
EXS=EX*EX 
EXY=£X*£Y 

Cl (K) sGMUJAC*  <FRT*tXS*EYS) 

C2  (K ) *GMUJAC*EXY*03 
C3(K)=GMUJAC*(£XS*FRT*EYS) 

C4(K)=GPRK*DJAC*(EXS*£YS> 

TC(K)=Q(JtK.4)*Rl-0,5*(U(K)**2*V(K)**2> 

CS1 (K)*GMUJAC*0Y 
CSS=-T3*CS1 <K)*V(K) 

CS2(K)=CS1 (K)*EX 
CS3  (K )  =CS1  (K)*tY 
CS4 (K) *-CSS* (U(K)*tX*V(K)*£Y) 

CS7  (K)*GPRK*0JAC*EY 

31  CONTINUE 

00  41  K=2»KM 

C5(K)s(TC(K*l)-TC(K-ln*0.S 
C6(K)=(  U(K»1)-  U(K-1 ) )*0.5 
C7(k>=(  V  (A* 1 ) -  V(K-1)>*0.5 
CS5(K)=CS2(K)*(U<K)*C7<K)-T3*Y(tc»*C6(K>> 
CS4(K)=CS3(K)*(U(K)*C5(K)*FRT»V(K)*C7(h  > ) 

41  CONTINUE 

CS()  )s(-3.*TC(l)  *4.*TC<2>-TC(3>  >*0.5 
CMl)sf(-3.*  U 1 1 )  *4»*  J(?>-  U(3)  )*0.5 
C7(l)*(-3.«  V)  1 ) *4. *  V<2)-  V ( 3) ) *0.5 
C5(KMAX)=(3.*TC(KM4X)-4.*TC  (Kh) *TC fH  M-l » ) *0.5 
C6 (KMAX) = (3.*  U(KMAX)-4.*  U  (AM) ♦  U <KH- 1 ) ) *0 .5 
C7  (KMAX)  =  (3.*  V(KMAXJ-4.*  V  (KM)»  V  t  K  *4- 1 )  )  *0.5 
00  32  K=2tKM 
KP=K*1 
KH=K— 1 

CSKPKP=CS1 (KP)*T3*Y(J»KP) 

C5KRKH=CS1 <KK)*T3*Y(J«KR) 

SP2=C1 (KP)*Cb(KP) ♦C2(XP)*C7(KP)-r5KPKP*XEX( J.K.2)*V(KP) 

SR2=C1 (KR)*Cb(KH) *C2 ( <R) *C7 (KR) -C5KRKR*XtX (J»K*2)*V(KR) 
FF(K«2>=SP2-SR2 

SP3=C2(KP)*Cb(KP>*C3(<P)*C7(KP>-r5KPKP*XtY(J.K,2)*V(KP) 

SR3=C2 (KP) *C6  (KR) »C3 ( KR) *C7 (KR) -C5KRKR*Xt Y ( J.K ,2) *V (KR) 

EF (K«3)=SP3-SR3 

SP4=C4  (KP) *C5 (KP) ♦ (Cl (KP)*U(KP)*r2<KP)*V(KP>  >*Cb(KP) 

1  ♦(C2(KP)*U<KP> ♦C3(KP)*V(KP>  >*C7(KP)-CS4(KP) 

SR4=C4  (KR)  *C5  (KR) ♦ (Cl  (KR)  *U  (KR)  *C.2  (KR)  *V  (KP)  )  *CM  KR> 

1  ♦ (C2(KR)*U(KH)*C3(RH)*V(KR) > *C7 (KR) -CS4 (KR) 

EF (K*4)=SP4-SR4 

T2S  <  CS3 (K ) *Cb (K)*CS2(K)*C7(K) )*?. 

T3=(2.*(CS3(K)*C7(k)-CS1 (K) *V (K) /V ( J.K) ) ) *2. 

T4* (CS7 (K)*C5(K)*(CS5(K)*C56(K)-T3*V ( K) **2/Y (J*K))*CSl (K))*2. 
EF  <K,2)*EF (K»2) *T2 
EF  (K *3)  *EF  (K  *3)  *T3 
EF(K.4)=EF(K*4)*T4 

32  CONTINUE 

00  33  K»2»KM 
00  33  N*2»4 

33  S(J.K.N)*S(J»K.N)*tF(K,N> 
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30  CONTINUE 
frFTlJNN 
END 
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Table  D4. 

Input  Data  Cards 

Card  No. 

Format 

Parameters 

1 

815 

JMAX,  KMAX, 

ITER,  IPRT,  IR1, 

2 

815 

JNM,  IGEOM, 

LIP,  KRES,  ITRAN 

3 

7F10.0, 

XMACH,  GAM, 

TM,  OMEGA, 

2F5.0 

CN,  CF,  CC, 

SMU,  SUMIMP 

IW2,  IFABD 
,  IVIS 


*  *  *  *  *  Stop  here  for  inviscid  flow  over  sphere-cone  and  follow  by  *  *  *  *  * 

last  Input  card,  Example  1 

4  5F10.0,  REY,  PRD,  PRT,  CVIS, 

315  TWA,  I TWA,  ITUR,  ITF 

it  *  *  *  *  Stop  here  for  viscous  (laminar)  flow  over  sphere-cone  and  *  *  *  * 

follow  by  last  Input  card.  Example  2 


Data  cards  1-3  are  always  needed  for  any  runs  of  Inviscid 
*  *  *  *  *  calculation.  Card  4  must  be  added  for  any  runs  of  viscous 

calculation 


*  *  *  *  * 


The  following  data  cards  are  needed  for  doing  arbitrary 
*  *  *  *  *  nosetip  shapes  with  cone  afterbody.  For  simplicity, 

examples  are  given  for  inviscid  flow  only.  For  viscous 
flows  the  same  additional  data  cards  must  be  included. 


*  *  *  *  *  IGE0M  =  read  in  XB>  YB,  XS,  YS  for  arbitrary 

nosetips  shape,  Example  3 

(4)l  8F10.5  XB,  YB,  XS,  YS  for  each  J 


*  *  *  *  * 


(4  +  JHAX)1 

Followed  by  the  last  input  card 


If  IGEOM  -  2,  read  in  Th(J)  and  DETT(J)  and  LIP  for  the 
*  *  *  *  *  nosetip  to  deform  from  sphere-cone  to  the  desired  nosetip 

shape,  Example  4 


*  *  *  *  * 
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I 

I 

k*  , 

I 

. 

J 


(4)2  8F10.5  TH(J) ,  J  -  1, JMAX 

(4  +  JMAX/8)2 

(5+JMAX/8)2  8F10.5  DETT(J) ,  J  -  1,  JMAX 

(5  +  (JMAX/6)x2)2 

Followed  by  the  last  input  card 


*  *  *  *  *  ^  IGEOM  =  3  or  4,  read  in  control  points  for  the  nosetip  *  *  *  *  * 

shape. 


<*>3.4 

7F10.0 

XI,  X2,  X3,  X4, 

X5, 

X6, 

X7 

<5)3,4 

7F10.0 

Yl,  Y2 ,  Y3,  Y4, 

Y5, 

Y6, 

Y7 

(6)3,4 

7F10.0 

Rl,  R2,  R3,  R4 

(8>3,4 

7F10.0 

SL1,  SL2 ,  SL3 

*  *  *  *  *  For  I^EOM  *  3  uniformly  distributed  TH(J)  in  the  nosetip  ***** 

portion,  then  the  last  input  card  follows.  Example  5. 

*  *  *  *  *  If  IGEOM  *  4,  read  in  TH(J)  in  the  nosetip  portion,  then  *  *  *  *  * 

the  last  input  card  follows.  Example  6. 

(9)4  8F10.0  TH(J) ,  J  -  1,  JMAX 


(9  +  JM^ 


*  *  *  *  * 


Note:  when  only  a  fraction  of  the  total  deformation  is  done  *  *  * 
for  this  run,  read  in  FACTB  before  the  last  input  card. 


n-A7 
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(9) 


(10  + 


3 

or 
JMAX 


8F10.0 


FACTB 


8  }4 


*  *  *  *  *  *or  a<*ding  points.'  OMEGA  =  0  and  IR1  =  1  after  reading  *  *  *  *  * 
cards  1,  2  and  3  the  4th  card  Is 


*  *  *  *  * 


4 

315 

JAA,  JIX,  KIM 

* 

When  more  rays  are  needed, 
the  THs  for  each 

JAA  0  and  KIM  =  0.  Read  in  ^  ^  ^ 
added  ray.  Example  7 

5 

8F10.0 

THAD(J) ,  J  -  1,  JAA 

Followed  by  reading  control  points  for  nosetlp  shape: 

6 

7F10.0 

XI,  X2 . .  X7 

7 

7F10.0 

Yl,  Y2,  . Y7 

8 

7F10.0 

Rl . .  R4 

9 

7F10.0 

SL1 y • « • • SL3 

Followed  by  the  last  Input  card. 

*  *  *  *  *  When  only  Increase  or  change  points  in  the  K  -  direction  ***** 

Then  JAA  =  0  and  KIM  t  0.  Example  8 


5’ 


F10.0  CF1 

Followed  by  the  last  input  card. 


Last  input  card  controls  the  output  of  flow  variables  at 
*****  each  K  line.  LP(K)  =  1  prints  out  the  flow  variables;  ***** 

LP(K)  *  0  skip  printing 


Last 


8011 


LP (K) ,  K  =  1,  KMAX 


When  starting  data  is  needed  for  afterbody  calculation,  let 
*****  IAFBD  *1.  A  data  card  for  JWRIT  must  be  read  in  after  the  *  *  *  *  * 

last  card.  Example  8. 


Final 


15 


JWRIT 


I 
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